Commit ac730c59 authored by Alain O' Miniussi's avatar Alain O' Miniussi
Browse files

Merge branch '184-remove-omeganew' into 'master'

Resolve "Remove OmegaNew"

Closes #184

See merge request DISC/fargOCA!103

Former-commit-id: c4d8f0e91e8c4f3b4f229f48098e551e96125001
parents 6e8b1548 3c75d4bf
......@@ -40,7 +40,6 @@ namespace fargOCA {
DR_ROTATING = 0x1, // possibly with omegaframe == 0 (could be FIXED)
DR_COROTATING = 0x2, // match the speed of the first planet
DR_GUIDING_CENTER = 0x4, // better approximation of co-rotating
DR_NEW_ROTATING = 0x8 // Moving from one speed to another.
};
inline
......@@ -139,8 +138,7 @@ namespace fargOCA {
struct Referential {
DiskReferentialType type = DR_NONE;
boost::optional<real> omega = boost::none;
boost::optional<real> newOmega = boost::none;
real omega = 0;
bool indirectForces = true;
bool trackPlanet = false;
TransportType transport = DT_NORMAL;
......@@ -149,7 +147,6 @@ namespace fargOCA {
bool rotating() const { return isType(DR_ROTATING); }
bool coRotating() const { return isType(DR_COROTATING); }
bool guidingCenter() const { return isType(DR_GUIDING_CENTER); }
bool newRotating() const { return isType(DR_NEW_ROTATING); }
bool fastTransport() const { return transport == DT_FAST; }
};
......
......@@ -154,7 +154,7 @@ main(int argc,
if (sys->restarted()) {
sys->restart("planet%d.dat", *restartOutput);
disk.grid().log() << "Planets restarted Omegaframe= " << *disk.physic().referential.omega
disk.grid().log() << "Planets restarted Omegaframe= " << disk.physic().referential.omega
<< "\t PhysicalTime = " << progress.physicalTime() << '\n';
} else {
if (disk.grid().master()) {
......
......@@ -96,9 +96,9 @@ namespace ADC {
/* centered-in-cell azimuthal velocity */
real vt_cent;
if (j < ns-1) {
vt_cent = 0.5*(vtheta[l]+vtheta[l+1]) + Rmed[i]**disk.physic().referential.omega;
vt_cent = 0.5*(vtheta[l]+vtheta[l+1]) + Rmed[i]*disk.physic().referential.omega;
} else {
vt_cent = 0.5*(vtheta[l]+vtheta[i*ns]) + Rmed[i]**disk.physic().referential.omega;
vt_cent = 0.5*(vtheta[l]+vtheta[i*ns]) + Rmed[i]*disk.physic().referential.omega;
}
total += Surf[i]*density[l]*Rmed[i]*vt_cent;
}
......@@ -143,7 +143,7 @@ namespace ADC {
/ (Rsup[i]-Rinf[i]));
/* centered-in-cell azimuthal velocity */
int nextl = j < ns-1 ? l+1 : i*ns;
real vt_cent = 0.5*(vtheta[l]+vtheta[nextl]) + Rmed[i]**disk.physic().referential.omega;
real vt_cent = 0.5*(vtheta[l]+vtheta[nextl]) + Rmed[i]*disk.physic().referential.omega;
total += Surf[i]*(0.5*density[l]*(vr_cent*vr_cent
+ vt_cent*vt_cent)
+ (energy - density[l]*pot[l]));
......@@ -196,7 +196,7 @@ namespace ADC {
real rot = ((Rmed[i+1]*vt[lip]-Rmed[im]*vt[lim]+Rmed[i+1]*vt[ljpip]-Rmed[im]*vt[limjp])
/ (Rmed[i]*2*(Rmed[i+1]-Rmed[im]))
+ 2* *disk.physic().referential.omega
+ 2* disk.physic().referential.omega
- (vr[ljp]-vr[ljm]+vr[ljpip]-vr[ljmip])*inv4dphi*InvRmed[i]);
vort[l] = rot/rho[l];
}
......
......@@ -115,7 +115,7 @@ namespace ADC {
real dx = Xplanet-xc;
real dy = Yplanet-yc;
real distance = sqrt(dx*dx+dy*dy);
real vtcell=0.5*(vtheta[l]+vtheta[ljp])+Rmed[i]**disk.physic().referential.omega;
real vtcell=0.5*(vtheta[l]+vtheta[ljp])+Rmed[i]*disk.physic().referential.omega;
real vrcell=0.5*(vrad[l]+vrad[lip]);
real vxcell=(vrcell*xc-vtcell*yc)/Rmed[i];
real vycell=(vrcell*yc+vtcell*xc)/Rmed[i];
......
......@@ -191,7 +191,7 @@ namespace ADC {
work in a frame non-centered on the star */
if (disk.physic().referential.coRotating()) {
real newOmega = sys.frequency(*corotMark, 0) / dt;
real domega = newOmega-*disk.physic().referential.omega;
real domega = newOmega-disk.physic().referential.omega;
CorrectVtheta (Vtheta, domega);
DiskPhysic::Referential ref = disk.physic().referential;
ref.omega = newOmega;
......@@ -199,7 +199,7 @@ namespace ADC {
initGrid.refreshAzimutalSpeed();
}
if (!disk.physic().referential.trackPlanet) {
sys.rotate(*disk.physic().referential.omega*dt);
sys.rotate(disk.physic().referential.omega*dt);
}
/* Now we update gas */
waveKiller(*Density, *Vrad, *Vtheta, dt);
......@@ -511,7 +511,7 @@ namespace ADC {
/* gravitational forces and curvature terms */
#pragma omp parallel
{
real omega = *disk.physic().referential.omega;
real omega = disk.physic().referential.omega;
#pragma omp for nowait
for (int i=1; i<nr; i++) {
for (int j=0; j<ns; j++) {
......
......@@ -78,7 +78,7 @@ namespace ADC {
real *tp = ThetaMomP ->field();
real *tm = ThetaMomM ->field();
auto const& Rmed = grid.rMed();
real omega = *disk.physic().referential.omega;
real omega = disk.physic().referential.omega;
#pragma omp parallel for
for (int i=0; i<nr; i++) {
for (int j=0; j<ns-1; j++) {
......@@ -467,7 +467,7 @@ namespace ADC {
} else if (i > 1 || (i > 0)) {
vr[l] = (rp[lim]+rm[l])/(rho[l]+rho[lim]+1e-20);
}
vt[l] = (tp[ljm]+tm[l])/(rho[l]+rho[ljm]+1e-15)/Rmed[i]-Rmed[i]**disk.physic().referential.omega;
vt[l] = (tp[ljm]+tm[l])/(rho[l]+rho[ljm]+1e-15)/Rmed[i]-Rmed[i]*disk.physic().referential.omega;
/* It was the angular momentum */
}
}
......
......@@ -89,7 +89,7 @@ namespace D3 {
vr[l+ni*ns] = vr[l+2*ni*ns];
vr[l] = 0.;
// Keplerian speed
vtheta[l] = (sqrt(G*1.0/Rmed[0])- *disk.physic().referential.omega*Rmed[0])*sin(phiMed[h]);
vtheta[l] = (sqrt(G*1.0/Rmed[0])- disk.physic().referential.omega*Rmed[0])*sin(phiMed[h]);
}
}
}
......@@ -113,7 +113,7 @@ namespace D3 {
vr[l+ni*ns] = -vr[l-ni*ns]; // Out of bounds
vr[l]=0;
// Check Bert code -> keplerian speed
vtheta[l] = (sqrt(G*1.0/Rmed[i])- *disk.physic().referential.omega*Rmed[i])*sin(phiMed[h]);
vtheta[l] = (sqrt(G*1.0/Rmed[i]) - disk.physic().referential.omega*Rmed[i])*sin(phiMed[h]);
}
}
}
......@@ -159,7 +159,7 @@ namespace D3 {
}
}
for (int h =0; h<ni; h++){
vtheta[ring(0,h)] = (sqrt(G*1/Rmed[0]) - *disk.physic().referential.omega*Rmed[0])*sin(phiMed[h]);
vtheta[ring(0,h)] = (sqrt(G*1/Rmed[0]) - disk.physic().referential.omega*Rmed[0])*sin(phiMed[h]);
}
}
......@@ -180,7 +180,7 @@ namespace D3 {
}
}
for (int h =0; h<ni; h++){
vtheta[ring(nr-1,h)] = (sqrt(G*1.0/Rmed[nr-1])- *disk.physic().referential.omega*Rmed[nr-1])*sin(phiMed[h]);
vtheta[ring(nr-1,h)] = (sqrt(G*1.0/Rmed[nr-1])- disk.physic().referential.omega*Rmed[nr-1])*sin(phiMed[h]);
}
}
}
......
......@@ -291,7 +291,7 @@ namespace D3 {
real gradp = (Press[l]-Press[lim])*2.0/(rho[l]+rho[lim])*invDiffRmed[i];
real gradphi = (Pot[l]-Pot[lim])*invDiffRmed[i];
real vt2 = std::pow((vtheta[l]+vtheta[ljp]+vtheta[l-ns*ni]+vtheta[ljp-ns*ni])/4
+ Rinf[i]**disk().physic().referential.omega*sinPhiMed[h],
+ Rinf[i]*disk().physic().referential.omega*sinPhiMed[h],
2);
real vp2 = std::pow(h < ni-1 // \todo recode loops for speed
? (vphi[l]+vphi[lhp]+vphi[lim]+vphi[lim+ns])/4
......@@ -326,7 +326,7 @@ namespace D3 {
real gradp = (Press[l]-Press[lhm])*2.0/(rho[l]+rho[lhm])*invDPhiMed[k];
real gradphi = (Pot[l]-Pot[lhm])*invDPhiMed[k];
real vt2 = vtheta[l]+vtheta[ljp]+vtheta[lhm]+vtheta[ljp-ns];
vt2 = vt2/4.0+Rmed[i]**disk().physic().referential.omega*sinPhi[h];
vt2 = vt2/4.0+Rmed[i]*disk().physic().referential.omega*sinPhi[h];
vt2 = vt2*vt2;
vphiint[l] = vphi[l]+dt*(-gradp-gradphi+vt2*cosPhi[h]/sinPhi[h]/Rmed[i]);
}
......@@ -521,21 +521,9 @@ namespace D3 {
if (!disk.physic().referential.trackPlanet) {
gas.system().advanceRK(dt);
}
if (disk.physic().referential.newRotating() && !system().restarted()){
real massTaper = progress().physicalTime()/(system().physic().massTaper.taper.final*2*PI);
massTaper = massTaper > 1 ? 1 : pow(sin(massTaper*PI/2.0),2.0);
DiskPhysic::Referential ref = disk.physic().referential;
real omega = *ref.omega;
real newOmega = *ref.newOmega;
real omegaNew = (newOmega-omega)*massTaper+omega;
real domega = omegaNew-omega;
CorrectVtheta(gas.velocity().theta(), domega);
ref.omega = omegaNew;
const_cast<DiskPhysic&>(disk.physic()).resetReferential(ref);
}
if (disk.physic().referential.coRotating()) {
DiskPhysic::Referential ref = disk.physic().referential;
real omega = *ref.omega;
real omega = ref.omega;
real newOmega = system().frequency(*corotMark, 0) / dt;
real domega = newOmega-omega;
CorrectVtheta(gas.velocity().theta(), domega);
......@@ -543,7 +531,7 @@ namespace D3 {
const_cast<DiskPhysic&>(disk.physic()).resetReferential(ref);
}
if(!disk.physic().referential.trackPlanet) {
gas.system().rotate(*disk.physic().referential.omega*dt);
gas.system().rotate(disk.physic().referential.omega*dt);
}
gas.applyBoundaryConditions(dt);
gas.computeSoundSpeed();
......
......@@ -431,7 +431,7 @@ namespace D3 {
auto const& Rmed = grid.rMed();
auto const& sinPhiMed = grid.sinPhiMed();
real omega = *disk().physic().referential.omega;
real omega = disk().physic().referential.omega;
#pragma omp parallel for
for (int i = 0; i < nr-1; i++) {
......@@ -508,7 +508,7 @@ namespace D3 {
if (j == 0)
ljm = h*ns+i*ni*ns+ns-1;
vr[l] = (rp[lim]+rm[l])/(rho[l]+rho[lim]);
vt[l] = (tp[ljm]+tm[l])/(rho[l]+rho[ljm])/(Rmed[i]*sinPhiMed[h])-Rmed[i]**disk().physic().referential.omega*sinPhiMed[h]; /* It was the angular momentum */
vt[l] = (tp[ljm]+tm[l])/(rho[l]+rho[ljm])/(Rmed[i]*sinPhiMed[h])-Rmed[i]*disk().physic().referential.omega*sinPhiMed[h]; /* It was the angular momentum */
vp[l] = (pp[lhm]+pm[l])/(rho[l]+rho[lhm])/Rmed[i];
}
}
......
......@@ -113,16 +113,12 @@ namespace fargOCA {
bool rotating = (t&DR_ROTATING) != 0;
bool coRotating = (t&DR_COROTATING) != 0;
bool guidingCenter = (t&DR_GUIDING_CENTER) != 0;
bool newRotating = (t&DR_NEW_ROTATING) != 0;
std::string str;
std::vector<std::string> mask;
if (guidingCenter) {
mask.push_back("GuidingCenter");
}
if (newRotating){
mask.push_back("NewRotating");
}
if (coRotating) {
mask.push_back("CoRotating");
}
......@@ -162,8 +158,6 @@ namespace fargOCA {
t = DR_COROTATING;
} else if (tk == "GUIDINGCENTER") {
t = DR_GUIDING_CENTER;
} else if (tk == "NEWROTATING") {
t = DR_NEW_ROTATING;
} else if (tk == "NONE") {
t = DR_NONE;
} else {
......
......@@ -265,7 +265,7 @@ namespace fargOCA {
real total = 0;
int nrMin = grid.managedRadiusMin(1);
int nrMax = grid.managedRadiusMax(1);
real omega = *disk().physic().referential.omega;
real omega = disk().physic().referential.omega;
#pragma omp parallel for reduction(+:total)
for (int i = nrMin; i < nrMax; i++) {
for (int h = 1; h <ni-1; h++) {
......@@ -704,7 +704,7 @@ namespace fargOCA {
real dy = Yplanet-yc;
real dz = Zplanet-zc;
real distance = sqrt(dx*dx+dy*dy+dz*dz);
real vtcell = 0.5*(vtheta[l]+vtheta[ljp])+Rmed[i]**disk().physic().referential.omega;
real vtcell = 0.5*(vtheta[l]+vtheta[ljp])+Rmed[i]*disk().physic().referential.omega;
real vrcell = 0.5*(vrad[l]+vrad[lip]);
real vpcell = 0.5*(vphi[l]+vphi[lhp]);
real vxcell = cos(theta[j])*(vrcell*sinPhiMed[h]+vpcell*cosPhiMed[h])-vtcell*sin(theta[j]);
......
......@@ -1119,9 +1119,8 @@ namespace fargOCA {
r.indirectForces = cfg.get<bool>("IndirectForces", true);
r.trackPlanet = cfg.get<bool>("TrackPlanet", false);
r.omega = cfg.get<real>("Omega", 0);
r.newOmega = cfg.get<real>("NewOmega", 0);
r.transport = cfg.get<TransportType>("Transport", DT_FAST);
if(r.trackPlanet && *r.omega == 0) {
if(r.trackPlanet && r.omega == 0) {
log << "You cannot use a fixed planet in a nonrotating frame\n"
<< "Either use Fixed=NO or use Fixed=YES and Frame ROTATING with "
<< "the *disk.physic().referential.omega corresponding to Keplerian velocity at planet position\n";
......
......@@ -429,7 +429,7 @@ namespace fargOCA {
<< ' ' << p.mass()
// 3D used to multiply those with RHO0*VOL0/XMSOL
<< ' ' << lostMass() << ' ' << lostMassZ() << ' ' << lostAccretion()
<< ' ' << time << ' ' << *(disk().physic().referential.omega)
<< ' ' << time << ' ' << disk().physic().referential.omega
<< '\n';
} else {
out << tag << '\t'
......
......@@ -43,7 +43,7 @@ namespace fargOCA {
return (KeplerianOmega(rad)*rad*sqrt(1 - (std::pow(disk.aspectRatio(rad),2.0)
* (1 + disk.physic().density.slope
- 2*disk.physic().flaringIndex) ))
- *disk.physic().referential.omega*rad);
- disk.physic().referential.omega*rad);
}
real
......@@ -98,7 +98,7 @@ namespace fargOCA {
G*(1/ri-1/r) / (r-ri));
}
real omega = *disk.physic().referential.omega;
real omega = disk.physic().referential.omega;
for (int i=1; i<grid.globalNbRadius(); i++) {
vt_int[i] = sqrt(vt_int[i]*grid.globalRadInf(i))-grid.globalRadInf(i)*omega;
}
......@@ -174,7 +174,7 @@ namespace fargOCA {
bool
TheoreticalProfiles2D::validAzimutalSpeed() const {
return std::abs(myLastOmega - *(disk().physic().referential.omega)) < 10*std::numeric_limits<real>::epsilon();
return std::abs(myLastOmega - disk().physic().referential.omega) < 10*std::numeric_limits<real>::epsilon();
}
void
......@@ -203,7 +203,7 @@ namespace fargOCA {
// Accessor will redirect unbalance request to the default array.
std::swap(myUnbalancedAzimutalSpeed, myAzimutalSpeed);
}
myLastOmega = *(disk().physic().referential.omega);
myLastOmega = disk().physic().referential.omega;
}
std::valarray<real> const&
......@@ -230,7 +230,7 @@ namespace fargOCA {
TheoreticalProfiles2D::TheoreticalProfiles2D(DiskDesc const& disk)
: TheoreticalProfiles(disk),
myLastOmega(*(disk.physic().referential.omega)),
myLastOmega(disk.physic().referential.omega),
myRadialSpeed(std::numeric_limits<real>::quiet_NaN(), disk.grid().nbRadius()),
myAzimutalSpeed(),
myUnbalancedAzimutalSpeed(),
......@@ -490,7 +490,7 @@ namespace fargOCA {
+ 2*flaringIndex
- (2*flaringIndex-1)*rg/Rmed[i]),
0.5)
- *disk().physic().referential.omega*Rmed[i]*sin(phiMed[h]));
- disk().physic().referential.omega*Rmed[i]*sin(phiMed[h]));
}
}
}
......
......@@ -34,10 +34,9 @@ main() {
DR_ROTATING,
DR_COROTATING,
DR_GUIDING_CENTER,
DR_NEW_ROTATING,
DR_GUIDING_CENTER|DR_NEW_ROTATING,
DR_GUIDING_CENTER|DR_NEW_ROTATING|DR_NONE,
DR_GUIDING_CENTER|DR_NEW_ROTATING|DR_ROTATING,
DR_GUIDING_CENTER,
DR_GUIDING_CENTER|DR_NONE,
DR_GUIDING_CENTER|DR_ROTATING,
};
bool passed = true;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment