r_mag added as a parameter for initializing mad B fields, and multizo…

…ne parfile updated

kharma/prob/seed_B.cpp

@@ -216,6 +216,7 @@ TaskStatus SeedBFieldType(MeshBlockData<Real> *rc, ParameterInput *pin, IndexDom
         // Require and load what we need if necessary
         Real A0 = pin->GetOrAddReal("b_field", "A0", 0.);
         Real min_A = pin->GetOrAddReal("b_field", "min_A", 0.2);
+        Real r_mag = pin->GetOrAddReal("b_field", "r_mag", 400.);
         // Init-specific loads
         Real a, rin, rmax, gam, kappa, rho_norm, arg1, n, rs, rb;
         Real tilt = 0; // Needs to be initialized
@@ -307,7 +308,7 @@ TaskStatus SeedBFieldType(MeshBlockData<Real> *rc, ParameterInput *pin, IndexDom
                     }
                 }
 
-                Real Aphi = seed_a<Seed>(Xmidplane, dxc, rho_av, rin, min_A, A0, arg1, rb);
+                Real Aphi = seed_a<Seed>(Xmidplane, dxc, rho_av, rin, min_A, A0, arg1, rb, r_mag);
 
                 if (tilt != 0.0) {
                     // This is *covariant* A_mu of an untilted disk

kharma/prob/seed_B.hpp

@@ -54,7 +54,7 @@ TaskStatus NormalizeBField(MeshData<Real> *md, ParameterInput *pin);
 enum BSeedType{constant, monopole, orszag_tang, orszag_tang_a, wave, shock_tube,
                 sane, mad, mad_quadrupole, r3s3, r5s5, gaussian, bz_monopole, vertical, r1s2, r1gizmo};
 
-#define SEEDA_ARGS GReal *x, const GReal *dxc, double rho, double rin, double min_A, double A0, double arg1, double rb
+#define SEEDA_ARGS GReal *x, const GReal *dxc, double rho, double rin, double min_A, double A0, double arg1, double rb, double r_mag
 
 // This will also act as the default implementation for unspecified types,
 // which should all be filled as B field by seed_b below.
@@ -81,7 +81,7 @@ template<>
 KOKKOS_INLINE_FUNCTION Real seed_a<BSeedType::mad>(SEEDA_ARGS)
 {
     return m::max(m::pow(x[1] / rin, 3) * m::pow(sin(x[2]), 3) *
-            m::exp(-x[1] / 400) * rho - min_A, 0.);
+            m::exp(-x[1] / r_mag) * rho - min_A, 0.);
 }
 
 // MAD, but turned into a quadrupole

pars/multizone/multizone.par

@@ -10,7 +10,7 @@ archive_parameters = true
 # Full mesh size, no refinement
 # Don't bother with xN boundaries for spherical coordinate systems
 # KHARMA will automatically place ~5 zones inside the EH
-nx1 = 96
+nx1 = 288
 nx2 = 64
 nx3 = 64
 
@@ -22,39 +22,62 @@ nx3 = 64
 
 <coordinates>
 base = ks
-transform = mks
-a = 0.0
-hslope = 0.3
-r_out = 512
+transform = eks
+a = 0.5
+r_out = 134217728
 r_in = 1
 
 <parthenon/time>
-tlim = 2000.0
+tlim = 10000000000
 
 <GRMHD>
-cfl = 0.9
-gamma = 1.666667
-reconstruction = weno5
+cfl = 0.45
+gamma = 1.6666666666666667
+
+<ismr>
+on = 1
+nlevels = 4
 
 <bondi>
-mdot = 1.0
-rs = 8.0
+mdot = 1.832747665572202
+rs = 316.22776601683796
+ur_frac = 1.0
+uphi = 0.0
+diffinit = true
+set_outer_bound = 0
+r_in_bondi = 2.0
 
 # Disable floors
 <floors>
-disable_floors = true
+disable_floors = false
+rho_min_geom = 1e-6
+u_min_geom = 1e-8
+bsq_over_u_max = 1000
+u_over_rho_max = 100
+bsq_over_rho_max = 100
+gamma_max = 5
+adjust_k = 0
+radius_dependent_gamma_max = 1
 
 # We'll be adding material, and that's okay
 <boundaries>
+inner_x1 = outflow
+outer_x1 = outflow
 check_inflow_outer_x1 = false
+check_inflow_inner_x1 = true
+
+<perturbation>
+jitter_above_rho = 0 # always apply perturbation when u_jitter > 0
+u_jitter = 0.1
 
 <b_field>
 solver = face_ct
 ct_scheme = sg07
 # Constant pure-vertical field
 type = r1s2
-A0 = 0.005
+A0 = 2e-07
 norm = false
+kill_on_large_divb = 0
 
 <debug>
 verbose = 1
@@ -64,24 +87,38 @@ extra_checks = 1
 <driver>
 type = multizone
 two_sync = true
+reconstruction = weno5
 
 <multizone>
 on = 1
-nzones = 2
+nzones = 8
+loc_tchar = 1
+ncycle_per_zone = 200
+one_trun = 1
+combine_out = false
+long_t_in = 10
+bflux_const = 1
+move_rin = 1
 
 <parthenon/output0>
 file_type = hdf5
-#dt = 5.0
-single_precision_output = false
-# Fields not present are silently ignored
-variables = prims, pflag, fflag, divB, cons.fB
+dn = -1
+single_precision_output = true
+variables = prims, cons.fB, fflag, pflag, divB
 ghost_zones = true
+analysis_output = 1
 
 <parthenon/output1>
 file_type = rst
-dn = 1000
+dn = 10000
 ghost_zones = true
 
 <parthenon/output2>
 file_type = hst
-dt = 0.1
+dn = 1000000
+
+<inverter>
+type = onedw
+
+<flux>
+consistent_face_b = 1