From 2d89f59215e81e7acd4f8062517894c0e8053d14 Mon Sep 17 00:00:00 2001 From: "yjiang@flatironinstitute.org" Date: Fri, 22 Nov 2024 09:36:11 -0500 Subject: [PATCH 1/5] fix the issue that default boundary conditions for specific intensities are not working for implicit radiation module --- src/bvals/cc/nr_radiation/bvals_rad.cpp | 249 ++++++++++++++++++++++++ src/bvals/cc/nr_radiation/bvals_rad.hpp | 7 + src/task_list/im_rad_task_list.cpp | 5 + src/task_list/im_rad_task_list.hpp | 2 + src/task_list/im_radit_task_list.cpp | 5 +- 5 files changed, 267 insertions(+), 1 deletion(-) diff --git a/src/bvals/cc/nr_radiation/bvals_rad.cpp b/src/bvals/cc/nr_radiation/bvals_rad.cpp index ea51de9373..0c227be951 100644 --- a/src/bvals/cc/nr_radiation/bvals_rad.cpp +++ b/src/bvals/cc/nr_radiation/bvals_rad.cpp @@ -17,10 +17,14 @@ #include "../../../coordinates/coordinates.hpp" #include "../../../globals.hpp" #include "../../../hydro/hydro.hpp" +#include "../../../field/field.hpp" #include "../../../mesh/mesh.hpp" #include "../../../nr_radiation/radiation.hpp" #include "bvals_rad.hpp" +#include // stringstream +#include // endl + //---------------------------------------------------------------------------------------- //! \class RadiationBoundaryFunctions @@ -398,3 +402,248 @@ void RadBoundaryVariable::PolarBoundarySingleAzimuthalBlock() { } return; } + + +//---------------------------------------------------------------------------------------- +//! \fn void RadBoundaryVariable::ApplyRadPhysicalBoundaries() +// \brief Apply physical boundaries for specific intensities only +void RadBoundaryVariable::ApplyRadPhysicalBoundaries(const Real time, const Real dt){ + + MeshBlock *pmb = pmy_block_; + Coordinates *pco = pmb->pcoord; + BoundaryValues *pbval = pmb->pbval; + int bis = pmb->is - NGHOST, bie = pmb->ie + NGHOST, + bjs = pmb->js, bje = pmb->je, + bks = pmb->ks, bke = pmb->ke; + + // Extend the transverse limits that correspond to periodic boundaries as they are + // updated: x1, then x2, then x3 + if (!pbval->apply_bndry_fn_[BoundaryFace::inner_x2] && pmb->block_size.nx2 > 1) + bjs = pmb->js - NGHOST; + if (!pbval->apply_bndry_fn_[BoundaryFace::outer_x2] && pmb->block_size.nx2 > 1) + bje = pmb->je + NGHOST; + if (!pbval->apply_bndry_fn_[BoundaryFace::inner_x3] && pmb->block_size.nx3 > 1) + bks = pmb->ks - NGHOST; + if (!pbval->apply_bndry_fn_[BoundaryFace::outer_x3] && pmb->block_size.nx3 > 1) + bke = pmb->ke + NGHOST; + + Hydro *ph = pmb->phydro; + + Field *pf = nullptr; + if (MAGNETIC_FIELDS_ENABLED) { + pf = pmb->pfield; + } + + NRRadiation *prad = pmb->pnrrad; + RadBoundaryVariable *pradbvar = &(prad->rad_bvar); + + + // Apply boundary function on inner-x1 + if (pbval->apply_bndry_fn_[BoundaryFace::inner_x1]) { + SetRadPhysicalFunctions(pmb, pco, time, dt, + pmb->is, pmb->ie, bjs, bje, bks, bke, NGHOST, + BoundaryFace::inner_x1, ph->w, pf->b, prad->ir); + } + + // Apply boundary function on outer-x1 + if (pbval->apply_bndry_fn_[BoundaryFace::outer_x1]) { + SetRadPhysicalFunctions(pmb, pco, time, dt, + pmb->is, pmb->ie, bjs, bje, bks, bke, NGHOST, + BoundaryFace::outer_x1, ph->w, pf->b, prad->ir); + + } + + if (pmb->block_size.nx2 > 1) { // 2D or 3D + // Apply boundary function on inner-x2 and update W,bcc (if not periodic) + if (pbval->apply_bndry_fn_[BoundaryFace::inner_x2]) { + SetRadPhysicalFunctions(pmb, pco, time, dt, + bis, bie, pmb->js, pmb->je, bks, bke, NGHOST, + BoundaryFace::inner_x2, ph->w, pf->b, prad->ir); + } + + if ((NR_RADIATION_ENABLED || IM_RADIATION_ENABLED) && + (pbval->block_bcs[BoundaryFace::inner_x2] != BoundaryFlag::block)) { + if (prad->rotate_theta == 1) { + pradbvar->RotateHPi_InnerX2(time, dt, bis, bie, pmb->js, bks, bke, NGHOST); + } + }// end radiation + + // Apply boundary function on outer-x2 and update W,bcc (if not periodic) + if (pbval->apply_bndry_fn_[BoundaryFace::outer_x2]) { + SetRadPhysicalFunctions(pmb, pco, time, dt, + bis, bie, pmb->js, pmb->je, bks, bke, NGHOST, + BoundaryFace::outer_x2, ph->w, pf->b, prad->ir); + } + } + + if (pmb->block_size.nx3 > 1) { // 3D + bjs = pmb->js - NGHOST; + bje = pmb->je + NGHOST; + + // Apply boundary function on inner-x3 and update W,bcc (if not periodic) + if (pbval->apply_bndry_fn_[BoundaryFace::inner_x3]) { + + SetRadPhysicalFunctions(pmb, pco, time, dt, + bis, bie, bjs, bje, pmb->ks, pmb->ke, NGHOST, + BoundaryFace::inner_x3, ph->w, pf->b, prad->ir); + } + + if ((NR_RADIATION_ENABLED || IM_RADIATION_ENABLED) && + (pbval->block_bcs[BoundaryFace::inner_x3] != BoundaryFlag::block)) { + if (prad->rotate_phi == 1) { + pradbvar->RotateHPi_InnerX3(time, dt, bis, bie, bjs, bje, pmb->ks, NGHOST); + } else if (prad->rotate_phi == 2) { + pradbvar->RotatePi_InnerX3(time, dt, bis, bie, bjs, bje, pmb->ks,NGHOST); + } + } + + // Apply boundary function on outer-x3 and update W,bcc (if not periodic) + if (pbval->apply_bndry_fn_[BoundaryFace::outer_x3]) { + SetRadPhysicalFunctions(pmb, pco, time, dt, + bis, bie, bjs, bje, pmb->ks, pmb->ke, NGHOST, + BoundaryFace::outer_x3, ph->w, pf->b, prad->ir); + } + if ((NR_RADIATION_ENABLED || IM_RADIATION_ENABLED) && + (pbval->block_bcs[BoundaryFace::outer_x3] != BoundaryFlag::block)) { + if (prad->rotate_phi == 1) { + pradbvar->RotateHPi_OuterX3(time, dt, bis, bie, bjs, bje, pmb->ke, NGHOST); + } else if (prad->rotate_phi == 2) { + pradbvar->RotatePi_OuterX3(time, dt, bis, bie, bjs, bje, pmb->ke, NGHOST); + } + } + } + return; +} + + + + +void RadBoundaryVariable::SetRadPhysicalFunctions( + MeshBlock *pmb, Coordinates *pco, Real time, Real dt, + int il, int iu, int jl, int ju, int kl, int ku, int ngh, + BoundaryFace face, AthenaArray &prim, + FaceField &b, AthenaArray &ir) { + + NRRadiation *prad = pmb->pnrrad; + RadBoundaryVariable *pradbvar = &(prad->rad_bvar); + BoundaryValues *pbval = pmb->pbval; + + + if (pbval->block_bcs[face] == BoundaryFlag::user) { // user-enrolled BCs + pmy_mesh_->RadBoundaryFunc_[face](pmb,pco,prad,prim,b, ir,time,dt, + il,iu,jl,ju,kl,ku,NGHOST); + } + // KGF: this is only to silence the compiler -Wswitch warnings about not handling the + // "undef" case when considering all possible BoundaryFace enumerator values. If "undef" + // is actually passed to this function, it will likely die before that ATHENA_ERROR() + // call, as it tries to access block_bcs[-1] + std::stringstream msg; + msg << "### FATAL ERROR in SetRadPhysicalFunctions" << std::endl + << "face = BoundaryFace::undef passed to this function" << std::endl; + + + switch (pbval->block_bcs[face]) { + case BoundaryFlag::user: // handled above, outside loop over BoundaryVariable objs + break; + case BoundaryFlag::reflect: + switch (face) { + case BoundaryFace::undef: + ATHENA_ERROR(msg); + case BoundaryFace::inner_x1: + pradbvar->ReflectInnerX1(time, dt, il, jl, ju, kl, ku, NGHOST); + break; + case BoundaryFace::outer_x1: + pradbvar->ReflectOuterX1(time, dt, iu, jl, ju, kl, ku, NGHOST); + break; + case BoundaryFace::inner_x2: + pradbvar->ReflectInnerX2(time, dt, il, iu, jl, kl, ku, NGHOST); + break; + case BoundaryFace::outer_x2: + pradbvar->ReflectOuterX2(time, dt, il, iu, ju, kl, ku, NGHOST); + break; + case BoundaryFace::inner_x3: + pradbvar->ReflectInnerX3(time, dt, il, iu, jl, ju, kl, NGHOST); + break; + case BoundaryFace::outer_x3: + pradbvar->ReflectOuterX3(time, dt, il, iu, jl, ju, ku, NGHOST); + break; + } + break; + case BoundaryFlag::outflow: + switch (face) { + case BoundaryFace::undef: + ATHENA_ERROR(msg); + case BoundaryFace::inner_x1: + pradbvar->OutflowInnerX1(time, dt, il, jl, ju, kl, ku, NGHOST); + break; + case BoundaryFace::outer_x1: + pradbvar->OutflowOuterX1(time, dt, iu, jl, ju, kl, ku, NGHOST); + break; + case BoundaryFace::inner_x2: + pradbvar->OutflowInnerX2(time, dt, il, iu, jl, kl, ku, NGHOST); + break; + case BoundaryFace::outer_x2: + pradbvar->OutflowOuterX2(time, dt, il, iu, ju, kl, ku, NGHOST); + break; + case BoundaryFace::inner_x3: + pradbvar->OutflowInnerX3(time, dt, il, iu, jl, ju, kl, NGHOST); + break; + case BoundaryFace::outer_x3: + pradbvar->OutflowOuterX3(time, dt, il, iu, jl, ju, ku, NGHOST); + break; + } + break; + case BoundaryFlag::vacuum: // special boundary condition type for radiation + switch (face) { + case BoundaryFace::undef: + ATHENA_ERROR(msg); + case BoundaryFace::inner_x1: + pradbvar->VacuumInnerX1(time, dt, il, jl, ju, kl, ku, NGHOST); + break; + case BoundaryFace::outer_x1: + pradbvar->VacuumOuterX1(time, dt, iu, jl, ju, kl, ku, NGHOST); + break; + case BoundaryFace::inner_x2: + pradbvar->VacuumInnerX2(time, dt, il, iu, jl, kl, ku, NGHOST); + break; + case BoundaryFace::outer_x2: + pradbvar->VacuumOuterX2(time, dt, il, iu, ju, kl, ku, NGHOST); + break; + case BoundaryFace::inner_x3: + pradbvar->VacuumInnerX3(time, dt, il, iu, jl, ju, kl, NGHOST); + break; + case BoundaryFace::outer_x3: + pradbvar->VacuumOuterX3(time, dt, il, iu, jl, ju, ku, NGHOST); + break; + } + break; + case BoundaryFlag::polar_wedge: + switch (face) { + case BoundaryFace::undef: + ATHENA_ERROR(msg); + case BoundaryFace::inner_x2: + pradbvar->PolarWedgeInnerX2(time, dt, il, iu, jl, kl, ku, NGHOST); + break; + case BoundaryFace::outer_x2: + pradbvar->PolarWedgeOuterX2(time, dt, il, iu, ju, kl, ku, NGHOST); + break; + default: + std::stringstream msg_polar; + msg_polar << "### FATAL ERROR in SetRadPhysicalBoundary" << std::endl + << "Attempting to call polar wedge boundary function on \n" + << "MeshBlock boundary other than inner x2 or outer x2" + << std::endl; + ATHENA_ERROR(msg_polar); + } + break; + default: + std::stringstream msg_flag; + msg_flag << "### FATAL ERROR in SetRadPhysicalBoundary" << std::endl + << "No BoundaryPhysics function associated with provided\n" + << "block_bcs[" << face << "] = BoundaryFlag::" + << GetBoundaryString(pbval->block_bcs[face]) << std::endl; + ATHENA_ERROR(msg); + break; + } // end switch (block_bcs[face]) +} + diff --git a/src/bvals/cc/nr_radiation/bvals_rad.hpp b/src/bvals/cc/nr_radiation/bvals_rad.hpp index 0d8d5bd27c..bd0f2a1666 100644 --- a/src/bvals/cc/nr_radiation/bvals_rad.hpp +++ b/src/bvals/cc/nr_radiation/bvals_rad.hpp @@ -38,6 +38,13 @@ class RadBoundaryVariable : public CellCenteredBoundaryVariable { void SetBoundaries() override; + void ApplyRadPhysicalBoundaries(const Real time, const Real dt); + void SetRadPhysicalFunctions( + MeshBlock *pmb, Coordinates *pco, Real time, Real dt, + int il, int iu, int jl, int ju, int kl, int ku, int ngh, + BoundaryFace face, AthenaArray &prim, + FaceField &b, AthenaArray &ir); + // function for shearing box void AddRadShearForInit(); void ShearQuantities(AthenaArray &shear_cc_, bool upper) override; diff --git a/src/task_list/im_rad_task_list.cpp b/src/task_list/im_rad_task_list.cpp index 63a14804cd..7812ca2f88 100644 --- a/src/task_list/im_rad_task_list.cpp +++ b/src/task_list/im_rad_task_list.cpp @@ -92,6 +92,11 @@ TaskStatus IMRadTaskList::PhysicalBoundary(MeshBlock *pmb) { return TaskStatus::success; } +TaskStatus IMRadTaskList::RadPhysicalBoundary(MeshBlock *pmb) { + pmb->pnrrad->rad_bvar.ApplyRadPhysicalBoundaries(time,dt); + return TaskStatus::success; +} + TaskStatus IMRadTaskList::ProlongateBoundary(MeshBlock *pmb) { pmb->pbval->ProlongateBoundaries(time, dt, pmb->pbval->bvars_main_int); diff --git a/src/task_list/im_rad_task_list.hpp b/src/task_list/im_rad_task_list.hpp index 04de7d8827..c77f0c582b 100644 --- a/src/task_list/im_rad_task_list.hpp +++ b/src/task_list/im_rad_task_list.hpp @@ -51,6 +51,8 @@ class IMRadTaskList { TaskListStatus DoAllAvailableTasks(MeshBlock *pmb, TaskStates &ts); void DoTaskListOneStage(Real wght); TaskStatus PhysicalBoundary(MeshBlock *pmb); + // set physical boundary for radiation only + TaskStatus RadPhysicalBoundary(MeshBlock *pmb); TaskStatus ProlongateBoundary(MeshBlock *pmb); protected: diff --git a/src/task_list/im_radit_task_list.cpp b/src/task_list/im_radit_task_list.cpp index 3eed88f56c..2649eda7ef 100644 --- a/src/task_list/im_radit_task_list.cpp +++ b/src/task_list/im_radit_task_list.cpp @@ -77,7 +77,7 @@ void IMRadITTaskList::AddTask(const TaskID& id, const TaskID& dep) { } else if (id == RAD_PHYS_BND) { task_list_[ntasks].TaskFunc= static_cast - (&IMRadITTaskList::PhysicalBoundary); + (&IMRadITTaskList::RadPhysicalBoundary); } else if (id == PRLN_RAD_BND) { task_list_[ntasks].TaskFunc= static_cast @@ -226,6 +226,9 @@ TaskStatus IMRadITTaskList::CheckResidual(MeshBlock *pmb) { } void IMRadITTaskList::StartupTaskList(MeshBlock *pmb) { + if(pmb->pnrrad->nfreq > 1) + pmb->pnrrad->UserFrequency(pmb->pnrrad); + pmb->pnrrad->rad_bvar.StartReceiving(BoundaryCommSubset::radiation); if (pmy_mesh->shear_periodic) { pmb->pnrrad->rad_bvar.StartReceivingShear(BoundaryCommSubset::radiation); From 9b8cfcf9954022d0579925fd80a7b43cc6ec15a4 Mon Sep 17 00:00:00 2001 From: "yjiang@flatironinstitute.org" Date: Fri, 22 Nov 2024 10:26:31 -0500 Subject: [PATCH 2/5] Add the option to allow two polar angles that are aligned with radial direction in spherical polar coordinate --- src/nr_radiation/angulargrid.cpp | 71 +++++++++++++++++++++++++++----- src/nr_radiation/radiation.cpp | 5 +++ src/nr_radiation/radiation.hpp | 1 + 3 files changed, 67 insertions(+), 10 deletions(-) diff --git a/src/nr_radiation/angulargrid.cpp b/src/nr_radiation/angulargrid.cpp index bf2070e6a6..5a9faf0370 100644 --- a/src/nr_radiation/angulargrid.cpp +++ b/src/nr_radiation/angulargrid.cpp @@ -41,6 +41,14 @@ void NRRadiation::AngularGrid(int angle_flag, int nmu) { AthenaArray pmat, pinv, wpf; AthenaArray plab, pl; + if(polar_angle){ + msg << "### FATAL ERROR in function [InitialAngle]" << std::endl + << "Polar angle cannot be used with: "<< angle_flag << "\n "; + throw std::runtime_error(msg.str().c_str()); + } + + + int n_ang = nang/noct; mu2tmp.NewAthenaArray(nmu); @@ -310,6 +318,11 @@ void NRRadiation::AngularGrid(int angle_flag, int nzeta, int npsi) { if (n2z > 1) ndim = 2; if (n3z > 1) ndim = 3; + Real coszeta_max = 1.0; + if(polar_angle){ + coszeta_max = ((Real)nang-2.0)/(Real)nang; + } + // separate ghost zones and active zones // so that they can be compatible with different angular scheme if (nzeta > 0) { @@ -323,10 +336,10 @@ void NRRadiation::AngularGrid(int angle_flag, int nzeta, int npsi) { int zs = 0; // ze = 2*nzeta - 1; - Real dcoszeta = 1.0/nzeta; - coszeta_f(zs) = 1.0; + Real dcoszeta = coszeta_max/nzeta; + coszeta_f(zs) = coszeta_max; for (int i=1; iGetOrAddInteger("radiation","npsi",0); angle_flag = pin->GetOrAddInteger("radiation","angle_flag",0); + polar_angle = pin->GetOrAddInteger("radiation","polar_angle",0); rotate_theta=pin->GetOrAddInteger("radiation","rotate_theta",0); rotate_phi=pin->GetOrAddInteger("radiation","rotate_phi",0); @@ -193,6 +194,10 @@ NRRadiation::NRRadiation(MeshBlock *pmb, ParameterInput *pin): nang = n_ang * noct; + // need to add two angles along the polar direction + if(polar_angle) + nang += 2; + // frequency grid //frequency grid covers -infty to infty, default nfreq=1, means gray // integrated over all frequency diff --git a/src/nr_radiation/radiation.hpp b/src/nr_radiation/radiation.hpp index 8b8fbc422a..a8e9d12621 100644 --- a/src/nr_radiation/radiation.hpp +++ b/src/nr_radiation/radiation.hpp @@ -77,6 +77,7 @@ class NRRadiation { int nang, noct, n_fre_ang; // n_fre_ang=nang*nfreq int angle_flag; + int polar_angle; // variables related to the angular space transport int nzeta, npsi; AthenaArray coszeta_v, zeta_v_full, zeta_f_full, dzeta_v, dzeta_f, From 9282f89eb55d17d8664969d14d43b39250a3054b Mon Sep 17 00:00:00 2001 From: Kyle Gerard Felker Date: Tue, 10 Dec 2024 20:11:42 -0600 Subject: [PATCH 3/5] Style errors --- src/bvals/cc/nr_radiation/bvals_rad.cpp | 15 +++++-------- src/nr_radiation/angulargrid.cpp | 28 ++++++++++--------------- 2 files changed, 16 insertions(+), 27 deletions(-) diff --git a/src/bvals/cc/nr_radiation/bvals_rad.cpp b/src/bvals/cc/nr_radiation/bvals_rad.cpp index 0c227be951..fae2950d2c 100644 --- a/src/bvals/cc/nr_radiation/bvals_rad.cpp +++ b/src/bvals/cc/nr_radiation/bvals_rad.cpp @@ -11,20 +11,19 @@ // C headers // C++ headers +#include // stringstream +#include // endl // Athena++ headers #include "../../../athena.hpp" #include "../../../coordinates/coordinates.hpp" +#include "../../../field/field.hpp" #include "../../../globals.hpp" #include "../../../hydro/hydro.hpp" -#include "../../../field/field.hpp" #include "../../../mesh/mesh.hpp" #include "../../../nr_radiation/radiation.hpp" #include "bvals_rad.hpp" -#include // stringstream -#include // endl - //---------------------------------------------------------------------------------------- //! \class RadiationBoundaryFunctions @@ -407,8 +406,7 @@ void RadBoundaryVariable::PolarBoundarySingleAzimuthalBlock() { //---------------------------------------------------------------------------------------- //! \fn void RadBoundaryVariable::ApplyRadPhysicalBoundaries() // \brief Apply physical boundaries for specific intensities only -void RadBoundaryVariable::ApplyRadPhysicalBoundaries(const Real time, const Real dt){ - +void RadBoundaryVariable::ApplyRadPhysicalBoundaries(const Real time, const Real dt) { MeshBlock *pmb = pmy_block_; Coordinates *pco = pmb->pcoord; BoundaryValues *pbval = pmb->pbval; @@ -450,7 +448,6 @@ void RadBoundaryVariable::ApplyRadPhysicalBoundaries(const Real time, const Real SetRadPhysicalFunctions(pmb, pco, time, dt, pmb->is, pmb->ie, bjs, bje, bks, bke, NGHOST, BoundaryFace::outer_x1, ph->w, pf->b, prad->ir); - } if (pmb->block_size.nx2 > 1) { // 2D or 3D @@ -466,7 +463,7 @@ void RadBoundaryVariable::ApplyRadPhysicalBoundaries(const Real time, const Real if (prad->rotate_theta == 1) { pradbvar->RotateHPi_InnerX2(time, dt, bis, bie, pmb->js, bks, bke, NGHOST); } - }// end radiation + } // end radiation // Apply boundary function on outer-x2 and update W,bcc (if not periodic) if (pbval->apply_bndry_fn_[BoundaryFace::outer_x2]) { @@ -482,7 +479,6 @@ void RadBoundaryVariable::ApplyRadPhysicalBoundaries(const Real time, const Real // Apply boundary function on inner-x3 and update W,bcc (if not periodic) if (pbval->apply_bndry_fn_[BoundaryFace::inner_x3]) { - SetRadPhysicalFunctions(pmb, pco, time, dt, bis, bie, bjs, bje, pmb->ks, pmb->ke, NGHOST, BoundaryFace::inner_x3, ph->w, pf->b, prad->ir); @@ -646,4 +642,3 @@ void RadBoundaryVariable::SetRadPhysicalFunctions( break; } // end switch (block_bcs[face]) } - diff --git a/src/nr_radiation/angulargrid.cpp b/src/nr_radiation/angulargrid.cpp index 5a9faf0370..e8b422cd9e 100644 --- a/src/nr_radiation/angulargrid.cpp +++ b/src/nr_radiation/angulargrid.cpp @@ -41,14 +41,12 @@ void NRRadiation::AngularGrid(int angle_flag, int nmu) { AthenaArray pmat, pinv, wpf; AthenaArray plab, pl; - if(polar_angle){ + if (polar_angle) { msg << "### FATAL ERROR in function [InitialAngle]" << std::endl << "Polar angle cannot be used with: "<< angle_flag << "\n "; throw std::runtime_error(msg.str().c_str()); } - - int n_ang = nang/noct; mu2tmp.NewAthenaArray(nmu); @@ -319,7 +317,7 @@ void NRRadiation::AngularGrid(int angle_flag, int nzeta, int npsi) { if (n3z > 1) ndim = 3; Real coszeta_max = 1.0; - if(polar_angle){ + if (polar_angle) { coszeta_max = ((Real)nang-2.0)/(Real)nang; } @@ -348,15 +346,13 @@ void NRRadiation::AngularGrid(int angle_flag, int nzeta, int npsi) { for (int i=0; i Date: Tue, 10 Dec 2024 20:13:07 -0600 Subject: [PATCH 4/5] Style error --- src/bvals/cc/nr_radiation/bvals_rad.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bvals/cc/nr_radiation/bvals_rad.cpp b/src/bvals/cc/nr_radiation/bvals_rad.cpp index fae2950d2c..1f319c9b47 100644 --- a/src/bvals/cc/nr_radiation/bvals_rad.cpp +++ b/src/bvals/cc/nr_radiation/bvals_rad.cpp @@ -11,8 +11,8 @@ // C headers // C++ headers -#include // stringstream #include // endl +#include // stringstream // Athena++ headers #include "../../../athena.hpp" From e850be62b3095d456d5d3ff5b1eb1ff413e913c3 Mon Sep 17 00:00:00 2001 From: "yjiang@flatironinstitute.org" Date: Wed, 11 Dec 2024 00:30:52 -0500 Subject: [PATCH 5/5] add radiation boundary during initialization for implicit_radiation --- src/mesh/mesh.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 5046fc79a3..06c2ad4479 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -1769,6 +1769,8 @@ void Mesh::Initialize(int res_flag, ParameterInput *pin) { } } pbval->ApplyPhysicalBoundaries(time, 0.0, pbval->bvars_main_int); + if(IM_RADIATION_ENABLED) + pmb->pnrrad->rad_bvar.ApplyRadPhysicalBoundaries(time,0.0); // Perform 4th order W(U) pmb->peos->ConservedToPrimitiveCellAverage(ph->u, ph->w1, pf->b, ph->w, pf->bcc, pmb->pcoord, @@ -1792,6 +1794,8 @@ void Mesh::Initialize(int res_flag, ParameterInput *pin) { } pbval->ApplyPhysicalBoundaries(time, 0.0, pbval->bvars_main_int); + if(IM_RADIATION_ENABLED) + pmb->pnrrad->rad_bvar.ApplyRadPhysicalBoundaries(time,0.0); } // for radiation, calculate opacity and moments if (NR_RADIATION_ENABLED || IM_RADIATION_ENABLED) {