From e83c23781eda97f3102b27aefaf402f88adbcef1 Mon Sep 17 00:00:00 2001 From: Vedant Dhruv Date: Wed, 6 Dec 2023 16:35:20 -0600 Subject: [PATCH] vertical b-field --- kharma/prob/chakrabarti_torus.cpp | 2 +- kharma/prob/seed_B.cpp | 72 +++++++++++++++++++++++--- kharma/prob/seed_B.hpp | 38 +++++++++++++- pars/tori_3d/chakrabarti.par | 85 +++++++++++++++++++++++++++++++ 4 files changed, 187 insertions(+), 10 deletions(-) create mode 100644 pars/tori_3d/chakrabarti.par diff --git a/kharma/prob/chakrabarti_torus.cpp b/kharma/prob/chakrabarti_torus.cpp index fe597685..1bbb75e5 100644 --- a/kharma/prob/chakrabarti_torus.cpp +++ b/kharma/prob/chakrabarti_torus.cpp @@ -73,7 +73,7 @@ TaskStatus InitializeChakrabartiTorus(std::shared_ptr>& 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 index 843d2678..8d198a13 100644 --- 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 *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 *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 *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("gamma"); + gam = pmb->packages.Get("GRMHD")->Param("gamma"); rho_norm = pmb->packages.Get("GRMHD")->Param("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("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 *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 *rc, ParameterInput *pin, IndexDom } } - Real Aphi = seed_a(Xmidplane, dxc, rho_av, rin, min_A, A0, arg1); + Real Aphi = seed_a(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 *md, ParameterInput *pin) status = SeedBFieldType(rc, pin); } else if (b_field_type == "vertical") { status = SeedBFieldType(rc, pin); + } else if (b_field_type == "vertical_chakrabarti") { + status = SeedBFieldType(rc, pin); } else if (b_field_type == "orszag_tang") { status = SeedBFieldType(rc, pin); } else if (b_field_type == "orszag_tang_a") { diff --git a/kharma/prob/seed_B.hpp b/kharma/prob/seed_B.hpp index dea13857..440b5601 100644 --- a/kharma/prob/seed_B.hpp +++ b/kharma/prob/seed_B.hpp @@ -52,9 +52,10 @@ TaskStatus NormalizeBField(MeshData *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(SEEDA_ARGS) + std::cos(x[2] + arg1)); } +template<> +KOKKOS_INLINE_FUNCTION Real seed_a(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 new file mode 100644 index 00000000..318f664d --- /dev/null +++ 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 + + +archive_parameters_timestamp = true +problem_id = torus +torus_type = chakrabarti + + +refinement = none +numlevel = 1 +nx1 = 256 +nx2 = 128 +nx3 = 128 + + +nx1 = 128 +nx2 = 128 +nx3 = 64 + + +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 + + +tlim = 10000.0 + + +cfl = 0.7 +gamma = 1.666667 +reconstruction = weno5 + + +type = imex +two_sync = true + + +rin = 20.0 +rmax = 41.0 + + +u_jitter = 0.1 + + +type = vertical_chakrabarti +beta_min = 100. + + +rho_min_geom = 1e-6 +u_min_geom = 1e-8 +bsq_over_rho_max = 100 +u_over_rho_max = 2 + + +verbose = 1 +extra_checks = 1 +flag_verbose = 0 + + +on = false +ne = 1.e-4 +Tp = 10 + + +file_type = hdf5 +dt = 5.0 +single_precision_output = true +variables = prims.rho, prims.u, prims.uvec, prims.B, jcon, fflag, pflag + + +file_type = rst +dt = 100.0 +ghost_zones = true + + +file_type = hst +dt = 0.1