diff --git a/kharma/b_cd/b_cd.cpp b/kharma/b_cd/b_cd.cpp index 730d273e..c426c8b9 100644 --- a/kharma/b_cd/b_cd.cpp +++ b/kharma/b_cd/b_cd.cpp @@ -1,25 +1,25 @@ -/* +/* * File: b_cd.cpp - * + * * BSD 3-Clause License - * + * * Copyright (c) 2020, AFD Group at UIUC * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, this * list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from * this software without specific prior written permission. - * + * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE @@ -246,7 +246,7 @@ TaskStatus PostStepDiagnostics(const SimTime& tm, MeshData *md) void FillOutput(MeshBlock *pmb, ParameterInput *pin) { - auto rc = pmb->meshblock_data.Get().get(); + auto rc = pmb->meshblock_data.Get("base").get(); IndexDomain domain = IndexDomain::interior; int is = pmb->cellbounds.is(domain), ie = pmb->cellbounds.ie(domain); int js = pmb->cellbounds.js(domain), je = pmb->cellbounds.je(domain); diff --git a/kharma/b_ct/b_ct.cpp b/kharma/b_ct/b_ct.cpp index 2e018a09..95fb1414 100644 --- a/kharma/b_ct/b_ct.cpp +++ b/kharma/b_ct/b_ct.cpp @@ -1,25 +1,25 @@ -/* +/* * File: b_ct.cpp - * + * * BSD 3-Clause License - * + * * Copyright (c) 2020, AFD Group at UIUC * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, this * list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from * this software without specific prior written permission. - * + * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE @@ -582,7 +582,8 @@ void B_CT::CalcDivB(MeshData *md, std::string divb_field_name) void B_CT::FillOutput(MeshBlock *pmb, ParameterInput *pin) { - auto rc = pmb->meshblock_data.Get(); + // This is called after the step, use "base" container + auto rc = pmb->meshblock_data.Get("base"); const int ndim = pmb->pmy_mesh->ndim; if (ndim < 2) return; diff --git a/kharma/b_ct/b_ct_boundaries.cpp b/kharma/b_ct/b_ct_boundaries.cpp index 242f657c..c2b44752 100644 --- a/kharma/b_ct/b_ct_boundaries.cpp +++ b/kharma/b_ct/b_ct_boundaries.cpp @@ -35,8 +35,10 @@ #include "decs.hpp" #include "domain.hpp" +#include "floors.hpp" #include "grmhd.hpp" #include "grmhd_functions.hpp" +#include "inverter.hpp" #include "kharma.hpp" void B_CT::ZeroBoundaryEMF(MeshBlockData *rc, IndexDomain domain, const VariablePack &emfpack, bool coarse) @@ -285,12 +287,21 @@ void B_CT::ReconnectBoundaryB3(MeshBlockData *rc, IndexDomain domain, cons const int bdir = KBoundaries::BoundaryDirection(bface); const auto bname = KBoundaries::BoundaryName(bface); + // Pull cell-centered values, as we need to update fluid primitives // TODO standardize on passing Packs or Datas... - auto B_U = rc->PackVariables(std::vector{"cons.B"}); - auto B_P = rc->PackVariables(std::vector{"prims.B"}); + PackIndexMap prims_map, cons_map; + auto P = rc->PackVariables({Metadata::GetUserFlag("Primitive"), Metadata::Cell}, prims_map); + auto U = rc->PackVariables(std::vector{Metadata::Conserved, Metadata::Cell}, cons_map); + const VarMap m_u(cons_map, true), m_p(prims_map, false); const auto& G = pmb->coords; + const Real gam = pmb->packages.Get("GRMHD")->Param("gamma"); + + const Floors::Prescription floors = pmb->packages.Get("Floors")->Param("prescription"); + // Don't be fooled, this function does *not* support/preserve EMHD values + const EMHD::EMHD_parameters& emhd_params = EMHD::GetEMHDParameters(pmb->packages); + // Subtract the average B3 as "reconnection" IndexRange3 b = KDomain::GetRange(rc, domain, F3, coarse); IndexRange3 bi = KDomain::GetRange(rc, IndexDomain::interior, F3, coarse); @@ -318,12 +329,22 @@ void B_CT::ReconnectBoundaryB3(MeshBlockData *rc, IndexDomain domain, cons ); member.team_barrier(); - // Update cell-centered conserved & primitive B3. Not worth a separate BlockUtoP call - parthenon::par_for_inner(member, b.ks, b.ke-1, + // Update cell-centered conserved & primitive B, and cell primitive fluid variables, in the zones we touched + parthenon::par_for_inner(member, b.ks, b.ke-1, // iterate over all *cells* k [&](const int& k) { - B_P(V3, k, jf, i) = (fpack(F3, 0, k, jf, i) / G.gdet(Loci::face3, jf, i) + P(m_p.B3, k, jf, i) = (fpack(F3, 0, k, jf, i) / G.gdet(Loci::face3, jf, i) + fpack(F3, 0, k + 1, jf, i) / G.gdet(Loci::face3, jf, i)) / 2; - B_U(V3, k, jf, i) = B_P(V3, k, jf, i) * G.gdet(Loci::center, jf, i); + U(m_u.B3, k, jf, i) = P(m_p.B3, k, jf, i) * G.gdet(Loci::center, jf, i); + + // Recover primitive GRMHD variables from our modified U + Inverter::u_to_p(G, U, m_u, gam, k, jf, i, P, m_p, Loci::center, + floors, 8, 1e-8); + // Floor them + int fflag = Floors::apply_geo_floors(G, P, m_p, gam, k, jf, i, floors, floors, Loci::center); + // Recalculate U on anything we floored + if (fflag) + GRMHD::p_to_u(G, P, m_p, gam, k, jf, i, U, m_u, Loci::center); + } ); } diff --git a/kharma/b_flux_ct/b_flux_ct.cpp b/kharma/b_flux_ct/b_flux_ct.cpp index e3813149..9f410a54 100644 --- a/kharma/b_flux_ct/b_flux_ct.cpp +++ b/kharma/b_flux_ct/b_flux_ct.cpp @@ -1,25 +1,25 @@ -/* +/* * File: b_flux_ct.cpp - * + * * BSD 3-Clause License - * + * * Copyright (c) 2020, AFD Group at UIUC * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, this * list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from * this software without specific prior written permission. - * + * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE @@ -410,7 +410,7 @@ void ZeroBoundaryFlux(MeshData *md, IndexDomain domain, bool coarse) // Compare this section with calculation of emf3 in FluxCT: // these changes ensure that boundary emfs emf3(i,js,k)=0, etc. for (auto &pmb : pmesh->block_list) { - auto& rc = pmb->meshblock_data.Get(); + auto& rc = pmb->meshblock_data.Get(md->StageName()); auto& B_F = rc->PackVariablesAndFluxes(std::vector{"cons.B"}); if (domain == IndexDomain::inner_x2 && @@ -495,7 +495,7 @@ void Bflux0(MeshData *md, IndexDomain domain, bool coarse) // Compare this section with calculation of emf3 in FluxCT: // these changes ensure that boundary emfs emf3(i,js,k)=0, etc. for (auto &pmb : pmesh->block_list) { - auto& rc = pmb->meshblock_data.Get(); + auto& rc = pmb->meshblock_data.Get(md->StageName()); auto& B_F = rc->PackVariablesAndFluxes(std::vector{"cons.B"}); // "Bflux0" prescription for keeping divB~=0 on zone corners of the interior & exterior X1 faces @@ -555,7 +555,7 @@ IndexRange ValidDivBX1(MeshBlock *pmb) { // All user, physical (not MPI/periodic) boundary conditions in X1 will generate divB on corners // intersecting the interior & exterior faces. Don't report these zones, as we expect it. - const IndexRange ibl = pmb->meshblock_data.Get()->GetBoundsI(IndexDomain::interior); + const IndexRange ibl = pmb->meshblock_data.Get("base")->GetBoundsI(IndexDomain::interior); bool avoid_inner = (!pmb->packages.Get("B_FluxCT")->Param("fix_flux_inner_x1") && pmb->boundary_flag[BoundaryFace::inner_x1] == BoundaryFlag::user); bool avoid_outer = (!pmb->packages.Get("B_FluxCT")->Param("fix_flux_outer_x1") && @@ -677,7 +677,7 @@ void CalcDivB(MeshData *md, std::string divb_field_name) } void FillOutput(MeshBlock *pmb, ParameterInput *pin) { - auto rc = pmb->meshblock_data.Get().get(); + auto rc = pmb->meshblock_data.Get("base").get(); const int ndim = pmb->pmy_mesh->ndim; if (ndim < 2) return; diff --git a/kharma/boundaries/boundaries.cpp b/kharma/boundaries/boundaries.cpp index 8dadef61..ac289b79 100644 --- a/kharma/boundaries/boundaries.cpp +++ b/kharma/boundaries/boundaries.cpp @@ -192,7 +192,7 @@ std::shared_ptr KBoundaries::Initialize(ParameterInput *pin, std: // introduce divergence to the first physical zone. bool clean_face_B = pin->GetOrAddBoolean("boundaries", "clean_face_B_" + bname, (btype == "outflow")); params.Add("clean_face_B_"+bname, clean_face_B); - // Forcibly reconnect field loops that get trapped around the polar boundary. Probably not needed anymore. + // Forcibly reconnect field loops that get trapped around the polar boundary bool reconnect_B3 = pin->GetOrAddBoolean("boundaries", "reconnect_B3_" + bname, false); params.Add("reconnect_B3_"+bname, reconnect_B3); @@ -626,7 +626,7 @@ TaskStatus KBoundaries::FixFlux(MeshData *md) const IndexRange3 b1 = KDomain::GetRange(md, IndexDomain::interior, -1, 1); for (auto &pmb : pmesh->block_list) { - auto &rc = pmb->meshblock_data.Get(); + auto &rc = pmb->meshblock_data.Get(md->StageName()); for (int i = 0; i < BOUNDARY_NFACES; i++) { BoundaryFace bface = (BoundaryFace)i; diff --git a/kharma/coord_output/coord_output.cpp b/kharma/coord_output/coord_output.cpp index f8dbbfbe..bd1eed56 100644 --- a/kharma/coord_output/coord_output.cpp +++ b/kharma/coord_output/coord_output.cpp @@ -1,25 +1,25 @@ -/* +/* * File: coord_output.cpp - * + * * BSD 3-Clause License - * + * * Copyright (c) 2020, AFD Group at UIUC * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, this * list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from * this software without specific prior written permission. - * + * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE @@ -70,7 +70,7 @@ std::shared_ptr CoordinateOutput::Initialize(ParameterInput *pin, pkg->AddField("coords.r", m0); pkg->AddField("coords.th", m0); pkg->AddField("coords.phi", m0); - + // Metric pkg->AddField("coords.gcon", m2); pkg->AddField("coords.gcov", m2); @@ -100,7 +100,7 @@ TaskStatus CoordinateOutput::BlockUserWorkBeforeOutput(MeshBlock *pmb, Parameter { auto& globals = pmb->packages.Get("Globals")->AllParams(); if (!globals.Get("in_loop")) { - auto rc = pmb->meshblock_data.Get(); + auto rc = pmb->meshblock_data.Get("base"); PackIndexMap geom_map; auto Geom = rc->PackVariables({Metadata::GetUserFlag("CoordinateOutput")}, geom_map); diff --git a/kharma/current/current.cpp b/kharma/current/current.cpp index b29e9604..cdd1ad6f 100644 --- a/kharma/current/current.cpp +++ b/kharma/current/current.cpp @@ -1,25 +1,25 @@ -/* +/* * File: current.cpp - * + * * BSD 3-Clause License - * + * * Copyright (c) 2020, AFD Group at UIUC * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, this * list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from * this software without specific prior written permission. - * + * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE @@ -118,7 +118,7 @@ void Current::FillOutput(MeshBlock *pmb, ParameterInput *pin) { // The "preserve" container will only exist after we've taken a step, // catch that situation - auto& rc1 = pmb->meshblock_data.Get(); + auto& rc1 = pmb->meshblock_data.Get("base"); auto rc0 = rc1; // Avoid writing rc0's type when initializing. Still light. try { // Get the state at beginning of the step diff --git a/kharma/driver/simple_step.cpp b/kharma/driver/simple_step.cpp index ada9e217..3883e1c8 100644 --- a/kharma/driver/simple_step.cpp +++ b/kharma/driver/simple_step.cpp @@ -1,25 +1,25 @@ -/* +/* * File: simple_step.cpp - * + * * BSD 3-Clause License - * + * * Copyright (c) 2020, AFD Group at UIUC * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, this * list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from * this software without specific prior written permission. - * + * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE @@ -52,7 +52,7 @@ TaskCollection KHARMADriver::MakeSimpleTaskCollection(BlockList_t &blocks, int s // Allocate the fluid states ("containers") we need for each block for (auto& pmb : blocks) { - auto &base = pmb->meshblock_data.Get(); + auto &base = pmb->meshblock_data.Get("base"); if (stage == 1) { pmb->meshblock_data.Add("dUdt", base); for (int i = 1; i < integrator->nstages; i++) @@ -68,7 +68,7 @@ TaskCollection KHARMADriver::MakeSimpleTaskCollection(BlockList_t &blocks, int s TaskRegion &single_tasklist_per_pack_region = tc.AddRegion(num_partitions); for (int i = 0; i < num_partitions; i++) { auto &tl = single_tasklist_per_pack_region[i]; - // Container names: + // Container names: // '_full_step_init' refers to the fluid state at the start of the full time step (Si in iharm3d) // '_sub_step_init' refers to the fluid state at the start of the sub step (Ss in iharm3d) // '_sub_step_final' refers to the fluid state at the end of the sub step (Sf in iharm3d) @@ -86,7 +86,7 @@ TaskCollection KHARMADriver::MakeSimpleTaskCollection(BlockList_t &blocks, int s // Any package modifications to the fluxes. e.g.: // 1. CT calculations for B field transport // 2. Zero fluxes through poles - // etc + // etc auto t_fix_flux = tl.AddTask(t_fluxes, Packages::FixFlux, md_sub_step_init.get()); // Apply the fluxes to calculate a change in cell-centered values "md_flux_src" diff --git a/kharma/flux/reconstruction.hpp b/kharma/flux/reconstruction.hpp index 0897abe4..ce4559b3 100644 --- a/kharma/flux/reconstruction.hpp +++ b/kharma/flux/reconstruction.hpp @@ -79,9 +79,9 @@ KOKKOS_FORCEINLINE_FUNCTION double Median(double a, double b, double c) template KOKKOS_FORCEINLINE_FUNCTION void reconstruct(RECONSTRUCT_ONE_ARGS) {} +// TODO better defaults? As-is, forgetting to implement these causes problems! template KOKKOS_FORCEINLINE_FUNCTION void reconstruct_left(RECONSTRUCT_ONE_LEFT_ARGS) {} - template KOKKOS_FORCEINLINE_FUNCTION void reconstruct_right(RECONSTRUCT_ONE_RIGHT_ARGS) {} @@ -260,7 +260,7 @@ KOKKOS_FORCEINLINE_FUNCTION void reconstruct(RECONSTRUCT_ONE Real dq = x4 - x3; dq = mc(x3 - x2, dq, 2.0); - const Real alpha_lin = 2.0 * alpha_r * alpha_l / (alpha_r + alpha_l); + const Real alpha_lin = clip(2.0 * alpha_r * alpha_l / (alpha_r + alpha_l), 0.0, 1.0); rout = alpha_lin * rout + (1.0 - alpha_lin) * (x3 + 0.5 * dq); lout = alpha_lin * lout + (1.0 - alpha_lin) * (x3 - 0.5 * dq); } @@ -493,6 +493,18 @@ KOKKOS_FORCEINLINE_FUNCTION void reconstruct(const Real &q_im2, cons } } } +template<> +KOKKOS_FORCEINLINE_FUNCTION void reconstruct_left(RECONSTRUCT_ONE_LEFT_ARGS) +{ + Real null; + reconstruct(x1, x2, x3, x4, x5, lout, null); +} +template<> +KOKKOS_FORCEINLINE_FUNCTION void reconstruct_right(RECONSTRUCT_ONE_RIGHT_ARGS) +{ + Real null; + reconstruct(x1, x2, x3, x4, x5, null, rout); +} // Row-wise implementations diff --git a/kharma/grmhd/grmhd.cpp b/kharma/grmhd/grmhd.cpp index 831070ec..81e5060d 100644 --- a/kharma/grmhd/grmhd.cpp +++ b/kharma/grmhd/grmhd.cpp @@ -297,7 +297,7 @@ Real EstimateTimestep(MeshData *md) // TODO maybe split normal, ISMR timesteps? Excised pole/recalculated ctop too? double min_ndt = std::numeric_limits::max(); for (auto &pmb : pmesh->block_list) { - auto rc = pmb->meshblock_data.Get().get(); + auto rc = pmb->meshblock_data.Get(md->StageName()).get(); // We only need this block-wise to check boundary flags for ISMR, could special-case that const bool polar_inner_x2 = pmb->boundary_flag[BoundaryFace::inner_x2] == BoundaryFlag::user; const bool polar_outer_x2 = pmb->boundary_flag[BoundaryFace::outer_x2] == BoundaryFlag::user; @@ -629,7 +629,7 @@ void UpdateAveragedCtop(MeshData *md) auto pmesh = md->GetMeshPointer(); auto& params = pmesh->packages.Get("Boundaries")->AllParams(); for (auto &pmb : pmesh->block_list) { - auto &rc = pmb->meshblock_data.Get(); + auto &rc = pmb->meshblock_data.Get(md->StageName()); for (int i = 0; i < BOUNDARY_NFACES; i++) { BoundaryFace bface = (BoundaryFace)i; auto bname = KBoundaries::BoundaryName(bface); diff --git a/kharma/implicit/implicit.cpp b/kharma/implicit/implicit.cpp index e3629f0c..2566e0ac 100644 --- a/kharma/implicit/implicit.cpp +++ b/kharma/implicit/implicit.cpp @@ -296,8 +296,10 @@ TaskStatus Implicit::Step(MeshData *md_full_step_init, MeshData *md_ #endif const auto& G = U_full_step_init_all.GetCoords(b); - // Solver performance diagnostics + // Solver performance diagnostics. Real &solve_fail = solve_fail_all(b, 0, k, j, i); + // Reset (locally and in View) at first iteration! + if (iter == 1) solve_fail = SolverStatusR::converged; // Perform the solve only if it hadn't failed in any of the previous iterations. if (solve_fail != SolverStatusR::fail) { diff --git a/kharma/prob/post_initialize.cpp b/kharma/prob/post_initialize.cpp index 3b1a7e5b..2c381e63 100644 --- a/kharma/prob/post_initialize.cpp +++ b/kharma/prob/post_initialize.cpp @@ -1,25 +1,25 @@ -/* +/* * File: post_initialize.cpp - * + * * BSD 3-Clause License - * + * * Copyright (c) 2020, AFD Group at UIUC * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, this * list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from * this software without specific prior written permission. - * + * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE @@ -97,7 +97,7 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart) // since seeding may be based on density if (pin->GetOrAddBoolean("blob", "add_blob", false)) { for (auto &pmb : pmesh->block_list) { - auto rc = pmb->meshblock_data.Get(); + auto rc = pmb->meshblock_data.Get("base"); // This inserts only in vicinity of some global r,th,phi InsertBlob(rc.get(), pin); } @@ -172,7 +172,7 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart) } // The e- initialization is called during problem initialization, but we want an option - // to force it -- for example, if restarting an ideal GRMHD run, + // to force it -- for example, if restarting an ideal GRMHD run, if (pkgs.count("Electrons") && pin->GetOrAddBoolean("electrons", "reinitialize", false)) { std::cout << "Reinitializing electron temperatures!" << std::endl; Electrons::MeshInitElectrons(md.get(), pin); diff --git a/kharma/prob/problem.cpp b/kharma/prob/problem.cpp index fd11e08e..3e7acc1a 100644 --- a/kharma/prob/problem.cpp +++ b/kharma/prob/problem.cpp @@ -1,25 +1,25 @@ -/* +/* * File: problem.cpp - * + * * BSD 3-Clause License - * + * * Copyright (c) 2020, AFD Group at UIUC * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, this * list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from * this software without specific prior written permission. - * + * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE @@ -69,7 +69,7 @@ using namespace parthenon; void KHARMA::ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { - auto rc = pmb->meshblock_data.Get(); + auto rc = pmb->meshblock_data.Get("base"); auto prob = pin->GetString("parthenon/job", "problem_id"); // Required parameter Flag("ProblemGenerator_"+prob); // Also just print this, it's important diff --git a/pars/tests/mhdmodes.par b/pars/tests/mhdmodes.par index 77e62928..4c142f34 100644 --- a/pars/tests/mhdmodes.par +++ b/pars/tests/mhdmodes.par @@ -53,6 +53,7 @@ integrator = rk2 dt_min = 0.0001 +# Courant number, where limit of 2o scheme is 1.0 regardless of dimensionality cfl = 0.9 gamma = 1.333333 # Whether to evolve these variables with an @@ -63,7 +64,7 @@ implicit = false # Local Lax-Friedrichs fluxes # `hlle` selects Harten, Lax, van Leer & Einfeldt # No production problems use HLLE as it has failure modes -# we do not test for. +# we do not catch in the code. type = llf # Lots of reconstruction options, see the wiki reconstruction = weno5 diff --git a/pars/tori_3d/mad.par b/pars/tori_3d/mad.par index 809217c9..701fa536 100644 --- a/pars/tori_3d/mad.par +++ b/pars/tori_3d/mad.par @@ -11,15 +11,15 @@ problem_id = torus archive_parameters = timestamp -nx1 = 288 -nx2 = 128 -nx3 = 128 +nx1 = 128 +nx2 = 64 +nx3 = 64 -nx1 = 144 +nx1 = 128 nx2 = 64 # Cannot break up phi with transmitting boundaries -nx3 = 128 +nx3 = 64 base = spherical_ks @@ -54,7 +54,9 @@ gamma = 1.666667 type = llf -reconstruction = weno5 +# WENO5 Adaptive Order "At Home" -- +# WENO5Z with fallback to linear recon in difficult regions +reconstruction = weno5_linear type = kharma @@ -69,6 +71,10 @@ outer_x2 = transmitting # This reduces the polar "wake," but shrinks the timestep, # and has caused instabilities. YMMV. #excise_polar_flux = true +# Reconnect phi component of B at the pole, as if it were one zone +# Dramatically reduces instabilities +reconnect_B3_inner_x2 = true +reconnect_B3_outer_x2 = true frame = drift diff --git a/pars/tori_3d/sane.par b/pars/tori_3d/sane.par index eccdf7a4..715dc52c 100644 --- a/pars/tori_3d/sane.par +++ b/pars/tori_3d/sane.par @@ -52,7 +52,9 @@ gamma = 1.666667 type = llf -reconstruction = weno5 +# WENO5 Adaptive Order "At Home" -- +# WENO5Z with fallback to linear recon in difficult regions +reconstruction = weno5_linear type = kharma @@ -67,6 +69,10 @@ outer_x2 = transmitting # This reduces the polar "wake," but shrinks the timestep, # and has caused instabilities. YMMV. #excise_polar_flux = true +# Reconnect phi component of B at the pole, as if it were one zone +# Dramatically reduces instabilities +reconnect_B3_inner_x2 = true +reconnect_B3_outer_x2 = true # Now with FOFC it might work to set "frame = normal",