Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Radiation boundary conditions are not applied in implicit radiation #603

Open
3 tasks done
apbailey opened this issue Jul 22, 2024 · 6 comments · May be fixed by #639
Open
3 tasks done

Radiation boundary conditions are not applied in implicit radiation #603

apbailey opened this issue Jul 22, 2024 · 6 comments · May be fixed by #639

Comments

@apbailey
Copy link
Collaborator

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
@yanfeij
Copy link
Contributor

yanfeij commented Jul 22, 2024 via email

@apbailey
Copy link
Collaborator Author

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

@felker
Copy link
Contributor

felker commented Nov 20, 2024

@yanfeij is there a plan to open a branch at least with this fix?

@yanfeij
Copy link
Contributor

yanfeij commented Nov 20, 2024 via email

@yanfeij
Copy link
Contributor

yanfeij commented Dec 10, 2024

@felker The fix is now in branch nrradfix (with one other new feature I added). Shall I make a pull request to merge it?

@felker
Copy link
Contributor

felker commented Dec 10, 2024

yes that'd be great!

@felker felker changed the title radiation boundaries not applied in implicit radiation Radiation boundary conditions are not applied in implicit radiation Dec 11, 2024
@felker felker linked a pull request Dec 11, 2024 that will close this issue
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
3 participants
@felker @yanfeij @apbailey