From 39e7f2cbd7c86905ac0eae998d04df10f308cb97 Mon Sep 17 00:00:00 2001 From: Ben Prather Date: Wed, 9 Oct 2024 16:17:59 -0600 Subject: [PATCH] Internally consistent boundary logic Move handling of GRMHD vars at boundaries to Inverter, try to make comments clearer --- kharma/boundaries/boundaries.cpp | 2 +- kharma/electrons/electrons.cpp | 30 ++++++++++++++++-------------- kharma/grmhd/grmhd.cpp | 26 +++++++++++--------------- kharma/inverter/inverter.cpp | 14 +++++++++++--- 4 files changed, 39 insertions(+), 33 deletions(-) diff --git a/kharma/boundaries/boundaries.cpp b/kharma/boundaries/boundaries.cpp index b1452e96..8dadef61 100644 --- a/kharma/boundaries/boundaries.cpp +++ b/kharma/boundaries/boundaries.cpp @@ -562,7 +562,7 @@ void KBoundaries::ApplyBoundary(std::shared_ptr> &rc, IndexD // There are two modes of operation here: if (sync_prims) { // 1. Exchange/prolongate/restrict PRIMITIVE variables: (ImEx driver) - // Primitive variables and conserved B field are marked FillGhost + // Primitive variables (except B field!) are marked FillGhost // Explicitly run UtoP on B field, then PtoU on everything // TODO there should be a set of B field wrappers that dispatch this auto pkgs = pmb->packages.AllPackages(); diff --git a/kharma/electrons/electrons.cpp b/kharma/electrons/electrons.cpp index 0c572d2b..6b36ca0b 100644 --- a/kharma/electrons/electrons.cpp +++ b/kharma/electrons/electrons.cpp @@ -1,25 +1,25 @@ -/* +/* * File: electrons.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 @@ -198,6 +198,8 @@ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr

BlockUtoP = Electrons::BlockUtoP; pkg->BoundaryUtoP = Electrons::BlockUtoP; + // If we wanted to apply the domian boundaries to primitive Electrons variables + //pkg->DomainBoundaryPtoU = Electrons::BlockPtoU; return pkg; } @@ -361,8 +363,8 @@ TaskStatus ApplyElectronHeating(MeshBlockData *rc_old, MeshBlockData // Heat different electron passives based on different dissipation fraction models // Expressions here closely adapted (read: stolen) from implementation in iharm3d // courtesy of Cesar Diaz, see https://github.com/AFD-Illinois/iharm3d - - // In all of these the electron entropy stored value is the entropy conserving solution + + // In all of these the electron entropy stored value is the entropy conserving solution // and then when updated it becomes the energy conserving solution if (m_p.K_CONSTANT >= 0) { const Real fel = fel_const; @@ -457,7 +459,7 @@ TaskStatus ApplyElectronHeating(MeshBlockData *rc_old, MeshBlockData const Real lx2= pmb->packages.Get("GRMHD")->Param("lx2"); const Real edot= pmb->packages.Get("GRMHD")->Param("drive_edot"); GridScalar alfven_speed = rc->Get("alfven_speed").data; - + int Nx1 = pmb->cellbounds.ncellsi(IndexDomain::interior); int Nx2 = pmb->cellbounds.ncellsj(IndexDomain::interior); Real *dv0 = (Real*) malloc(sizeof(Real)*Nx1*Nx2); @@ -492,9 +494,9 @@ TaskStatus ApplyElectronHeating(MeshBlockData *rc_old, MeshBlockData dv0[i*Nx1+j] -= mean_velocity0; dv1[i*Nx1+j] -= mean_velocity1; } - } + } - Real Bhalf = 0; Real A = 0; Real init_e = 0; + Real Bhalf = 0; Real A = 0; Real init_e = 0; Kokkos::Sum Bhalf_reducer(Bhalf); Kokkos::Sum A_reducer(A); Kokkos::Sum init_e_reducer(init_e); pmb->par_reduce("forced_mhd_normal_kick_normalization_Bhalf", mykb.s, mykb.e, myjb.s, myjb.e, myib.s, myib.e, KOKKOS_LAMBDA(const int k, const int j, const int i, Real &local_result) { @@ -578,7 +580,7 @@ void ApplyFloors(MeshBlockData *mbd, IndexDomain domain) } else { ktot_max = floors.ktot_max; } - + if (P(m_p.KTOT, k, j, i) > ktot_max) { fflag(0, k, j, i) = Floors::FFlag::KTOT | (int) fflag(0, k, j, i); P(m_p.KTOT, k, j, i) = ktot_max; diff --git a/kharma/grmhd/grmhd.cpp b/kharma/grmhd/grmhd.cpp index dce07d80..5d119ece 100644 --- a/kharma/grmhd/grmhd.cpp +++ b/kharma/grmhd/grmhd.cpp @@ -1,25 +1,25 @@ -/* +/* * File: grmhd.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 @@ -190,16 +190,12 @@ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr

BlockUtoP // Taken care of by separate "Inverter" package since it's hard to do - - // On physical boundaries, even if we've sync'd both, respect the application to primitive variables - pkg->DomainBoundaryPtoU = Flux::BlockPtoUMHD; - - // AMR-related pkg->CheckRefinementBlock = GRMHD::CheckRefinement; pkg->EstimateTimestepMesh = GRMHD::EstimateTimestep; pkg->PostStepDiagnosticsMesh = GRMHD::PostStepDiagnostics; + // NOTE: PtoU and UtoP for the five fluid variables are taken care of by the "Inverter" + // package, which registers the relevant callbacks + // This is so that the whole package can be disabled when running implicitly or w/EMHD // List (vector) of HistoryOutputVars that will all be enrolled as output variables parthenon::HstVar_list hst_vars = {}; diff --git a/kharma/inverter/inverter.cpp b/kharma/inverter/inverter.cpp index 423c31f9..bdd3022d 100644 --- a/kharma/inverter/inverter.cpp +++ b/kharma/inverter/inverter.cpp @@ -107,9 +107,17 @@ std::shared_ptr Inverter::Initialize(ParameterInput *pin, std::sh m = Metadata({Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::OneCopy, Metadata::Overridable}); pkg->AddField("fflag", m); - // We exist basically to do this - pkg->BlockUtoP = Inverter::BlockUtoP; - pkg->BoundaryUtoP = Inverter::BlockUtoP; + // This package may be loaded even when evolving implicitly, e.g. for FOFC + // Only register our callbacks if they're needed for explicit evolution or a guess + if (!pin->GetBoolean("GRMHD", "implicit") || pin->GetBoolean("emhd", "ideal_guess")) { + // We exist basically to do this + pkg->BlockUtoP = Inverter::BlockUtoP; + // We want to run U->P on most boundaries when we're synchronizing conserved variables + pkg->BoundaryUtoP = Inverter::BlockUtoP; + // However, we apply domain boundaries to primitives. + // Registering this additional function conveys that to the callers in `Packages` and `Boundaries` + pkg->DomainBoundaryPtoU = Flux::BlockPtoUMHD; + } pkg->PostStepDiagnosticsMesh = Inverter::PostStepDiagnostics;