Cleaning

src/eos/adiabatic_glmmhd.hpp

@@ -93,9 +93,9 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
 
     // Let's apply floors explicitly, i.e., by default floor will be disabled (<=0)
     // and the code will fail if a negative density is encountered.
-    PARTHENON_REQUIRE(u_d > 0.0 || density_floor_ > 0.0,
-                      "Got negative density. Consider enabling first-order flux "
-                      "correction or setting a reasonable density floor.");
+    //PARTHENON_REQUIRE(u_d > 0.0 || density_floor_ > 0.0,
+    //                  "Got negative density. Consider enabling first-order flux "
+    //                  "correction or setting a reasonable density floor.");
 
     // apply density floor, without changing momentum or energy
     u_d = (u_d > density_floor_) ? u_d : density_floor_;
@@ -134,9 +134,9 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
 
     // Let's apply floors explicitly, i.e., by default floor will be disabled (<=0)
     // and the code will fail if a negative pressure is encountered.
-    PARTHENON_REQUIRE(w_p > 0.0 || pressure_floor_ > 0.0 || e_floor_ > 0.0,
-                      "Got negative pressure. Consider enabling first-order flux "
-                      "correction or setting a reasonable pressure or temperature floor.");
+    //PARTHENON_REQUIRE(w_p > 0.0 || pressure_floor_ > 0.0 || e_floor_ > 0.0,
+    //                  "Got negative pressure. Consider enabling first-order flux "
+    //                  "correction or setting a reasonable pressure or temperature floor.");
 
     // Pressure floor (if present) takes precedence over temperature floor
     if ((pressure_floor_ > 0.0) && (w_p < pressure_floor_)) {

src/hydro/srcterms/gravitational_field.hpp

@@ -36,7 +36,18 @@ void GravitationalFieldSrcTerm(parthenon::MeshData<parthenon::Real> *md,
   IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior);
   IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior);
   IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior);
-
+
+  /*
+  std::cout << "Values of gravitationnal field" << std::endl;
+  std::cout 
+    << gravitationalField.g_from_r(1e-4) << ", " << gravitationalField.g_from_r(3.16e-4)
+    << ", " << gravitationalField.g_from_r(1e-3) << ", " << gravitationalField.g_from_r(3.16e-3)
+    << ", " << gravitationalField.g_from_r(1e-2) << ", " << gravitationalField.g_from_r(3.16e-2) 
+    << ", " << gravitationalField.g_from_r(1e-1) << ", " << gravitationalField.g_from_r(3.16e-1)
+    << ", " << gravitationalField.g_from_r(1.0)  << ", " << gravitationalField.g_from_r(3.16) 
+    << ", " << gravitationalField.g_from_r(10.0) << std::endl;
+  */
+
   parthenon::par_for(
       DEFAULT_LOOP_PATTERN, "GravitationalFieldSrcTerm", parthenon::DevExecSpace(), 0,
       cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,

src/pgen/cloud.cpp

@@ -143,6 +143,7 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin) {
 //  \brief Problem Generator for the cloud in wind setup
 
 void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
+
   auto hydro_pkg = pmb->packages.Get("Hydro");
   auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
   auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
@@ -156,9 +157,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
     PARTHENON_FAIL("Requested to initialize magnetic fields by `cloud/plasma_beta > 0`, "
                    "but `hydro/fluid` is not supporting MHD.");
   }
-
+  
   auto steepness = pin->GetOrAddReal("problem/cloud", "cloud_steepness", 10);
-
+  
   // initialize conserved variables
   auto &mbd = pmb->meshblock_data.Get();
   auto &u_dev = mbd->Get("cons").data;
@@ -191,13 +192,13 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
         u(IM2, k, j, i) = mom;
         // Can use rhoe_wind here as simulation is setup in pressure equil.
         u(IEN, k, j, i) = rhoe_wind + 0.5 * mom * mom / rho;
-
+        
         if (mhd_enabled) {
           u(IB1, k, j, i) = Bx;
           u(IB2, k, j, i) = By;
           u(IEN, k, j, i) += 0.5 * (Bx * Bx + By * By);
         }
-
+        
         // Init passive scalars
         for (auto n = nhydro; n < nhydro + nscalars; n++) {
           if (rad <= r_cloud) {
@@ -240,13 +241,14 @@ void InflowWindX2(std::shared_ptr<MeshBlockData<Real>> &mbd, bool coarse) {
 }
 
 parthenon::AmrTag ProblemCheckRefinementBlock(MeshBlockData<Real> *mbd) {
+
   auto pmb = mbd->GetBlockPointer();
   auto w = mbd->Get("prim").data;
-
+  
   IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
   IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
   IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
-
+  
   auto hydro_pkg = pmb->packages.Get("Hydro");
   const auto nhydro = hydro_pkg->Param<int>("nhydro");