radiation boundaries not applied in implicit radiation

[apbailey]

Prerequisite checklist

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

• [X] Did you check that the issue hasn't already been reported?
• [X] Did you check the documentation in the Wiki for an answer?
• [X] 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

[yanfeij]

This is known and has been fixed but not committed to master yet. You can use user defined boundary. Or I can send you a version I have.

On Mon, Jul 22, 2024 at 6:41 PM Avery Bailey @.***> wrote:

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 https://github.com/user-attachments/files/16340332/pgen.txt input.txt https://github.com/user-attachments/files/16340333/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

— Reply to this email directly, view it on GitHub https://github.com/PrincetonUniversity/athena/issues/603, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACCWEJ2JFJ7I5XILAML2VJTZNWDCHAVCNFSM6AAAAABLJGEDOOVHI2DSMVQWIX3LMV43ASLTON2WKOZSGQZDGOBYGUZTEMY . You are receiving this because you are subscribed to this thread.Message ID: @.***>

-- Best,

Yan-Fei Jiang Associate Research Scientist, Flatiron Institute, 160 Fifth Avenue New York, NY 10010

[apbailey]

sounds good, just wanted to make sure it was known, I will just use user-defined until it gets committed to master