Shoot donuts into space

parthenon-hpc-lab · Aug 15, 2023 · bc0d804 · bc0d804
1 parent ea838bf
commit bc0d804
Show file tree

Hide file tree

Showing 2 changed files with 18 additions and 10 deletions.
diff --git a/src/pgen/cluster/magnetic_tower.cpp b/src/pgen/cluster/magnetic_tower.cpp
@@ -133,9 +133,9 @@ void MagneticTower::AddSrcTerm(parthenon::Real field_to_add, parthenon::Real mas
                               0.5 * (b_x * b_x + b_y * b_y + b_z * b_z);
 
         // Add density
-        const Real cell_delta_rho =
-            mt.DensityFromSimCart(coords.Xc<1>(i), coords.Xc<2>(j), coords.Xc<3>(k));
-        cons(IDN, k, j, i) += cell_delta_rho;
+        // const Real cell_delta_rho =
+        // mt.DensityFromSimCart(coords.Xc<1>(i), coords.Xc<2>(j), coords.Xc<3>(k));
+        // cons(IDN, k, j, i) += cell_delta_rho;
       });
 }
 

diff --git a/src/pgen/cluster/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp
@@ -50,11 +50,15 @@ class MagneticTowerObj {
   PotentialInJetCyl(const parthenon::Real r, const parthenon::Real h,
                     parthenon::Real &a_r, parthenon::Real &a_theta,
                     parthenon::Real &a_h) const __attribute__((always_inline)) {
-    const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2) - pow(h / l_scale_, 2));
+    const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2));
     // Compute the potential in jet cylindrical coordinates
     a_r = 0.0;
-    a_theta = field_ * l_scale_ * (r / l_scale_) * exp_r2_h2;
-    a_h = field_ * l_scale_ * alpha_ / 2.0 * exp_r2_h2;
+    a_theta = 0.0;
+    if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) {
+      a_h = field_ * l_scale_ * alpha_ / 2.0 * exp_r2_h2;
+    } else {
+      a_h = 0.0;
+    }
   }
 
   // Compute Magnetic Potential in simulation Cartesian coordinates
@@ -81,11 +85,15 @@ class MagneticTowerObj {
                 parthenon::Real &b_theta, parthenon::Real &b_h) const
       __attribute__((always_inline)) {
 
-    const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2) - pow(h / l_scale_, 2));
+    const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2));
     // Compute the field in jet cylindrical coordinates
-    b_r = field_ * 2 * (h / l_scale_) * (r / l_scale_) * exp_r2_h2;
-    b_theta = field_ * alpha_ * (r / l_scale_) * exp_r2_h2;
-    b_h = field_ * 2 * (1 - pow(r / l_scale_, 2)) * exp_r2_h2;
+    b_r = 0.0;
+    if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) {
+      b_theta = 2.0 * field_ * r / l_scale_ * exp_r2_h2;
+    } else {
+      b_theta = 0.0;
+    }
+    b_h = 0.0;
   }
 
   // Compute Magnetic field in Simulation Cartesian coordinates