diff --git a/src/bvals/cc/nr_radiation/bvals_rad.cpp b/src/bvals/cc/nr_radiation/bvals_rad.cpp index ea51de937..1f319c9b4 100644 --- a/src/bvals/cc/nr_radiation/bvals_rad.cpp +++ b/src/bvals/cc/nr_radiation/bvals_rad.cpp @@ -11,10 +11,13 @@ // C headers // C++ headers +#include // endl +#include // stringstream // Athena++ headers #include "../../../athena.hpp" #include "../../../coordinates/coordinates.hpp" +#include "../../../field/field.hpp" #include "../../../globals.hpp" #include "../../../hydro/hydro.hpp" #include "../../../mesh/mesh.hpp" @@ -398,3 +401,244 @@ 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 0d8d5bd27..bd0f2a166 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/mesh/mesh.cpp b/src/mesh/mesh.cpp index 5046fc79a..06c2ad447 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) { diff --git a/src/nr_radiation/angulargrid.cpp b/src/nr_radiation/angulargrid.cpp index bf2070e6a..e8b422cd9 100644 --- a/src/nr_radiation/angulargrid.cpp +++ b/src/nr_radiation/angulargrid.cpp @@ -41,6 +41,12 @@ 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 +316,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 +334,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 8b8fbc422..a8e9d1262 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, diff --git a/src/task_list/im_rad_task_list.cpp b/src/task_list/im_rad_task_list.cpp index 63a14804c..7812ca2f8 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 04de7d882..c77f0c582 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 3eed88f56..2649eda7e 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);