vertical b-field

AFD-Illinois · Dec 6, 2023 · e83c237 · e83c237
1 parent fd5ba95
commit e83c237
Show file tree

Hide file tree

Showing 4 changed files with 187 additions and 10 deletions.
diff --git a/kharma/prob/chakrabarti_torus.cpp b/kharma/prob/chakrabarti_torus.cpp
@@ -73,7 +73,7 @@ TaskStatus InitializeChakrabartiTorus(std::shared_ptr<MeshBlockData<Real>>& rc,
 
     // Thermodynamic quantities at rin, rmax
     const Real lnh_in             = lnh_calc(a, rin, rin, 1.0, cc, nn);
-    const Real lnh_peak            = lnh_calc(a, rin, rmax, 1.0, cc, nn) - lnh_in;
+    const Real lnh_peak           = lnh_calc(a, rin, rmax, 1.0, cc, nn) - lnh_in;
     const Real pgas_over_rho_peak = gm1/gam * (m::exp(lnh_peak) - 1.0);
     const Real rho_peak           = m::pow(pgas_over_rho_peak, 1.0 / gm1) / rho_max;
 

diff --git a/kharma/prob/seed_B.cpp b/kharma/prob/seed_B.cpp
@@ -40,6 +40,7 @@
 #include "coordinate_utils.hpp"
 #include "domain.hpp"
 #include "fm_torus.hpp"
+#include "chakrabarti_torus.hpp"
 #include "grmhd_functions.hpp"
 
 using namespace parthenon;
@@ -91,6 +92,17 @@ TaskStatus SeedBFieldType(MeshBlockData<Real> *rc, ParameterInput *pin, IndexDom
     std::string b_field_type = pin->GetString("b_field", "type");
     auto prob = pin->GetString("parthenon/job", "problem_id");
     bool is_torus = (prob == "torus");
+    // What kind?
+    std::string torus_type;
+    if (is_torus) {
+        torus_type = pin->GetString("parthenon/job", "torus_type");
+    }
+    bool is_fm, is_chakrabarti = false;
+    if (torus_type == "fishbone_moncrief") {
+        is_fm = true;
+    } else if (torus_type == "chakrabarti") {
+        is_chakrabarti = true;
+    }
 
     // Indices
     // TODO handle filling faces with domain < entire more gracefully
@@ -192,6 +204,11 @@ TaskStatus SeedBFieldType(MeshBlockData<Real> *rc, ParameterInput *pin, IndexDom
         // Init-specific loads
         Real a, rin, rmax, gam, kappa, rho_norm, arg1;
         Real tilt = 0; // Needs to be initialized
+        // Chakrabarti-specific variables
+        const Real rho_max = pin->GetOrAddReal("torus", "rho_max", 1.0);
+        Real gm1, lnh_in, lnh_peak, pgas_over_rho_peak, rho_peak;
+        GReal cc, nn; 
+        Real potential_rho_pow, potential_falloff, potential_r_pow;
         switch (Seed) {
         case BSeedType::sane:
         case BSeedType::mad:
@@ -200,14 +217,36 @@ TaskStatus SeedBFieldType(MeshBlockData<Real> *rc, ParameterInput *pin, IndexDom
         case BSeedType::r5s5:
         case BSeedType::gaussian:
             // Torus parameters
-            rin = pin->GetReal("torus", "rin");
-            rmax = pin->GetReal("torus", "rmax");
+            rin   = pin->GetReal("torus", "rin");
+            rmax  = pin->GetReal("torus", "rmax");
             kappa = pin->GetReal("torus", "kappa");
-            tilt = pin->GetReal("torus", "tilt") / 180. * M_PI;
+            tilt  = pin->GetReal("torus", "tilt") / 180. * M_PI;
             // Other things we need only for torus evaluation
-            gam = pmb->packages.Get("GRMHD")->Param<Real>("gamma");
+            gam      = pmb->packages.Get("GRMHD")->Param<Real>("gamma");
             rho_norm = pmb->packages.Get("GRMHD")->Param<Real>("rho_norm");
+            a        = G.coords.get_a();
+            break;
+        case BSeedType::vertical_chakrabarti:
+            // A separate case for the vertical field initialized for the Chakrabarti torus
+            // Fluid parameters
+            gam = pmb->packages.Get("GRMHD")->Param<Real>("gamma");
+            gm1 = gam - 1.;
+            // Field init parameters
+            potential_rho_pow = pin->GetOrAddReal("b_field", "potential_rho_pow", 1.0);
+            potential_falloff = pin->GetOrAddReal("b_field", "potential_falloff", 0.0);
+            potential_r_pow   = pin->GetOrAddReal("b_field", "potential_r_pow", 0.0);
+            // Torus parameters
+            rin  = pin->GetReal("torus", "rin");
+            rmax = pin->GetReal("torus", "rmax");
+            tilt = pin->GetReal("torus", "tilt") / 180. * M_PI;
+            // Spacetime geometry parameters
             a = G.coords.get_a();
+            // Now compute relevant parameters
+            cn_calc(a, rin, rmax, &cc, &nn);
+            lnh_in             = lnh_calc(a, rin, rin, 1.0, cc, nn);
+            lnh_peak           = lnh_calc(a, rin, rmax, 1.0, cc, nn) - lnh_in;
+            pgas_over_rho_peak = gm1/gam * (m::exp(lnh_peak) - 1.0);
+            rho_peak           = m::pow(pgas_over_rho_peak, 1.0 / gm1) / rho_max;
             break;
         case BSeedType::orszag_tang_a:
             A0 = pin->GetReal("orszag_tang", "tscale");
@@ -233,15 +272,31 @@ TaskStatus SeedBFieldType(MeshBlockData<Real> *rc, ParameterInput *pin, IndexDom
                 rotate_polar(Xembed, tilt, Xmidplane);
                 const GReal r = Xmidplane[1], th = Xmidplane[2];
 
+                // Trigonometric values
+                const GReal sth = sin(th);
+                const GReal cth = cos(th);
+
                 // In case we need zone sizes
                 const GReal dxc[GR_DIM] = {0., G.Dxc<1>(i), G.Dxc<2>(j), G.Dxc<3>(k)};
 
                 // This is written under the assumption re-computed rho is more accurate than a bunch
                 // of averaging in a meaningful way.  Just use the average if not.
                 Real rho_av;
+                bool in_torus = false;
                 if (is_torus) {
-                    // Find rho at corner directly for torii
-                    rho_av = fm_torus_rho(a, rin, rmax, gam, kappa, r, th) / rho_norm;
+                    if (is_fm) {
+                        // Find rho at corner directly for torii
+                        rho_av = fm_torus_rho(a, rin, rmax, gam, kappa, r, th) / rho_norm;
+                    }
+                    else if (is_chakrabarti){
+                        // Find rho
+                        const Real lnh = lnh_calc(a, rin, r, sth, cc, nn);
+                        if (lnh >= 0.0) {
+                            in_torus = true;
+                            Real pg_over_rho = gm1 / gam * (m::exp(lnh) - 1.0);
+                            rho_av = m::pow(pg_over_rho, 1. / gm1) / rho_peak;
+                        }
+                    }
                 } else {
                     // Use averages for anything else
                     // This loop runs over every corner. Centers do not exist before the first
@@ -265,7 +320,8 @@ TaskStatus SeedBFieldType(MeshBlockData<Real> *rc, ParameterInput *pin, IndexDom
                     }
                 }
 
-                Real Aphi = seed_a<Seed>(Xmidplane, dxc, rho_av, rin, min_A, A0, arg1);
+                Real Aphi = seed_a<Seed>(Xmidplane, dxc, rho_av, rin, min_A, A0, arg1, in_torus, rho_max,\
+                                        potential_rho_pow, potential_falloff, potential_r_pow);
 
                 if (tilt != 0.0) {
                     // This is *covariant* A_mu of an untilted disk
@@ -388,6 +444,8 @@ TaskStatus SeedBField(MeshData<Real> *md, ParameterInput *pin)
             status = SeedBFieldType<BSeedType::bz_monopole>(rc, pin);
         } else if (b_field_type == "vertical") {
             status = SeedBFieldType<BSeedType::vertical>(rc, pin);
+        } else if (b_field_type == "vertical_chakrabarti") {
+            status = SeedBFieldType<BSeedType::vertical_chakrabarti>(rc, pin);
         } else if (b_field_type == "orszag_tang") {
             status = SeedBFieldType<BSeedType::orszag_tang>(rc, pin);
         } else if (b_field_type == "orszag_tang_a") {

diff --git a/kharma/prob/seed_B.hpp b/kharma/prob/seed_B.hpp
@@ -52,9 +52,10 @@ TaskStatus NormalizeBField(MeshData<Real> *md, ParameterInput *pin);
 
 // Internal representation of the field initialization preference, used for templating
 enum BSeedType{constant, monopole, monopole_cube, orszag_tang, orszag_tang_a, wave, shock_tube,
-                sane, mad, mad_quadrupole, r3s3, r5s5, gaussian, bz_monopole, vertical};
+                sane, mad, mad_quadrupole, r3s3, r5s5, gaussian, bz_monopole, vertical, vertical_chakrabarti};
 
-#define SEEDA_ARGS GReal *x, const GReal *dxc, double rho, double rin, double min_A, double A0, double arg1
+#define SEEDA_ARGS GReal *x, const GReal *dxc, double rho, double rin, double min_A, double A0, double arg1, bool in_torus,\
+                            double rho_max, double potential_rho_pow, double potential_falloff, double potential_r_pow
 
 // This will also act as the default implementation for unspecified types,
 // which should all be filled as B field by seed_b below.
@@ -130,6 +131,39 @@ KOKKOS_INLINE_FUNCTION Real seed_a<BSeedType::orszag_tang_a>(SEEDA_ARGS)
                         + std::cos(x[2] + arg1));
 }
 
+template<>
+KOKKOS_INLINE_FUNCTION Real seed_a<BSeedType::vertical_chakrabarti>(SEEDA_ARGS)
+{
+    Real r      = x[1];
+    Real th     = x[2];
+    Real sth    = sin(x[2]);
+    Real svarth = std::fabs(sin(th));
+
+    Real cyl_radius   = r * svarth;
+    Real rcyl_in      = rin;
+    Real rcyl_falloff = potential_falloff;
+
+    Real Aphi = m::pow(cyl_radius / rcyl_in, potential_r_pow);
+    if (potential_falloff != 0) {
+        Aphi *= m::exp(-cyl_radius / rcyl_falloff);
+    }
+    Real Aphi_offset = exp(-rcyl_in/rcyl_falloff);
+    if (cyl_radius < rcyl_in) {
+      Aphi = 0.0;
+    } else {
+      Aphi -= Aphi_offset;
+    }
+    if (potential_rho_pow != 0.0) {
+        if (in_torus) {
+            Aphi *= m::pow(rho / rho_max, potential_rho_pow);
+        } else {
+            Aphi = 0.0;
+        }
+    }
+
+    return Aphi;
+}
+
 #undef SEEDA_ARGS
 #define SEEDB_ARGS GReal *x, GReal gdet, double k1, double k2, double k3, double phase, \
                     double amp_B1, double amp_B2, double amp_B3, \

diff --git a/pars/tori_3d/chakrabarti.par b/pars/tori_3d/chakrabarti.par
@@ -0,0 +1,85 @@
+# Chakrabarti torus (based on AthenaK implementation)
+# courtesy G. N. Wong.
+# b_field/type must be 'vertical_chakrabarti' (since there's a preexisting 'vertical' field type)
+# setup details
+
+<parthenon/job>
+archive_parameters_timestamp = true
+problem_id = torus
+torus_type = chakrabarti
+
+<parthenon/mesh>
+refinement = none
+numlevel = 1
+nx1 = 256
+nx2 = 128
+nx3 = 128
+
+<parthenon/meshblock>
+nx1 = 128
+nx2 = 128
+nx3 = 64
+
+<coordinates>
+base       = spherical_ks
+transform  = fmks
+r_out      = 1000
+a          = 0.9375
+hslope     = 0.3
+mks_smooth = 0.5
+poly_xt    = 0.82
+poly_alpha = 14.0
+
+<parthenon/time>
+tlim = 10000.0
+
+<GRMHD>
+cfl            = 0.7
+gamma          = 1.666667
+reconstruction = weno5
+
+<driver>
+type     = imex
+two_sync = true
+
+<torus>
+rin  = 20.0
+rmax = 41.0
+
+<perturbation>
+u_jitter = 0.1
+
+<b_field>
+type     = vertical_chakrabarti
+beta_min = 100.
+
+<floors>
+rho_min_geom     = 1e-6
+u_min_geom       = 1e-8
+bsq_over_rho_max = 100
+u_over_rho_max   = 2
+
+<debug>
+verbose      = 1
+extra_checks = 1
+flag_verbose = 0
+
+<wind>
+on = false
+ne = 1.e-4
+Tp = 10
+
+<parthenon/output0>
+file_type = hdf5
+dt = 5.0
+single_precision_output = true
+variables = prims.rho, prims.u, prims.uvec, prims.B, jcon, fflag, pflag
+
+<parthenon/output1>
+file_type = rst
+dt = 100.0
+ghost_zones = true
+
+<parthenon/output2>
+file_type = hst
+dt = 0.1