Prerequisite checklist

Place an X in between the brackets on each line as you complete these checks:

• Did you check that the issue hasn't already been reported?
• Did you check the documentation in the Wiki for an answer?
• Are you running the latest version of Athena++?

Summary of issue
It seems that with the implicit radiation module, physical radiation boundary conditions are never actually applied when bc=reflecting is used for example.

Steps to reproduce
pgen.txt
input.txt
I have attached a mwe problem generator and input file as txt files to reproduce this issue. The setup is static, homogenous 1D cartesian profile in radiative equilibrium, (rho=p=T=Er=intensity=opacity=1, v=0). When bc=reflecting, the equilibrium is not maintained, the reflecting BC functions are never called, the intensity in the ghost cells is zero and the gas cools. However, I have taken the code blocks for reflecting BCs and copied them into user-defined BCs so that if bc=user, a user defined reflecting BC is used -- this setup preserves equilibrium exactly. The equilibrium can also be maintained exactly with bc=reflecting, if one uses -nr_radiation instead of -implicit_radiation. While the ApplyPhysicalBoundaries function in the implicit loop im_rad_task_list.cpp appears to be called

TaskStatus IMRadTaskList::PhysicalBoundary(MeshBlock *pmb) {
  pmb->pnrrad->rad_bvar.var_cc = &(pmb->pnrrad->ir);
  pmb->pbval->ApplyPhysicalBoundaries(time, dt, pmb->pbval->bvars_main_int);
  return TaskStatus::success;
}

the RadBoundaryVariable::ReflectInnerX1 in src/bvals/cc/nr_radiation/reflect_rad.cpp is never called suggesting that the rad bcs are never actually enrolled in bvars_main_int. And indeed in radiation.cpp, push_back of the radiation bc is done to bvars, not bvars_main_int if using implicit radiation:

  pmb->pbval->bvars.push_back(&rad_bvar);
  // enroll radiation boundary value object
  if (NR_RADIATION_ENABLED) {
    pmb->pbval->bvars_main_int.push_back(&rad_bvar);
  }

Additional comments and/or proposed solution
There are two ways I have found to address this problem. The first works perfectly but it calls bvars_main_int in the implicit loop when maybe this is meant to be reserved for only the main integration loop so I do not know if it aligns with Yanfei/Athena intended design. The second works fairly well but doesn't fix to same level of precision.

Fix 1:
The first is to have the implicit iteration continue to apply physical boundaries to bvars_main_int but now enroll the rad bcs into bvars_main_int to mimic the -nr_radiation case. On my test problem, this fix gives identical results between bc=reflecting and bc=user. Here is the git diff

diff --git a/src/mesh/meshblock.cpp b/src/mesh/meshblock.cpp
index 3e3d0e41..13421f02 100644
--- a/src/mesh/meshblock.cpp
+++ b/src/mesh/meshblock.cpp
@@ -249,7 +249,7 @@ MeshBlock::MeshBlock(int igid, int ilid, LogicalLocation iloc, RegionSize input_
        //radiation constructor needs the parameter nfre_ang
     pnrrad = new NRRadiation(this, pin);
   }
-  if (NR_RADIATION_ENABLED) {
+  if (NR_RADIATION_ENABLED || IM_RADIATION_ENABLED) {
     pbval->AdvanceCounterPhysID(RadBoundaryVariable::max_phys_id);
   }
 
@@ -448,7 +448,7 @@ MeshBlock::MeshBlock(int igid, int ilid, Mesh *pm, ParameterInput *pin,
        //radiation constructor needs the parameter nfre_ang
     pnrrad = new NRRadiation(this, pin);
   }
-  if (NR_RADIATION_ENABLED) {
+  if (NR_RADIATION_ENABLED || IM_RADIATION_ENABLED) {
     pbval->AdvanceCounterPhysID(RadBoundaryVariable::max_phys_id);
   }
 
diff --git a/src/nr_radiation/radiation.cpp b/src/nr_radiation/radiation.cpp
index f0cf0f2b..ec27044e 100644
--- a/src/nr_radiation/radiation.cpp
+++ b/src/nr_radiation/radiation.cpp
@@ -292,7 +292,7 @@ NRRadiation::NRRadiation(MeshBlock *pmb, ParameterInput *pin):
   rad_bvar.bvar_index = pmb->pbval->bvars.size();
   pmb->pbval->bvars.push_back(&rad_bvar);
   // enroll radiation boundary value object
-  if (NR_RADIATION_ENABLED) {
+  if (NR_RADIATION_ENABLED || IM_RADIATION_ENABLED) {
     pmb->pbval->bvars_main_int.push_back(&rad_bvar);
   }

Fix 2
Instead of ApplyPhysicalBoundaries operating on bvars_main_int, have it operate on bvars where the radiation boundary is currently enrolled. On my test problem, this does not give identical results between bc=reflecting and bc=user. When bc=reflecting there are small ~1e-10 deviations from equilibrium, as compared with 0.0 deviation for bc=user. Perhaps there is another spot in the integration where ApplyPhysicalBoundaries is not operating on the radiation boundary values but I have not been able to find anything. Here is that git diff

diff --git a/src/task_list/im_rad_task_list.cpp b/src/task_list/im_rad_task_list.cpp
index 63a14804..9163f5f1 100644
--- a/src/task_list/im_rad_task_list.cpp
+++ b/src/task_list/im_rad_task_list.cpp
@@ -88,7 +88,7 @@ void IMRadTaskList::DoTaskListOneStage(Real wght) {
 
 TaskStatus IMRadTaskList::PhysicalBoundary(MeshBlock *pmb) {
   pmb->pnrrad->rad_bvar.var_cc = &(pmb->pnrrad->ir);
-  pmb->pbval->ApplyPhysicalBoundaries(time, dt, pmb->pbval->bvars_main_int);
+  pmb->pbval->ApplyPhysicalBoundaries(time, dt, pmb->pbval->bvars);
   return TaskStatus::success;
 }

Version info

• Athena++ version: 24.0 release
• Compiler and version: gcc