Skip to content

Commit

Permalink
Merge branch 'development' into sdc_clip
Browse files Browse the repository at this point in the history
  • Loading branch information
@zingale
zingale authored Oct 12, 2023
2 parents 1877ce2 + f8446cd commit 9328a6d
Show file tree
Hide file tree
Showing 19 changed files with 357 additions and 865 deletions.
28 changes: 26 additions & 2 deletions Docs/source/reactions.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,8 +57,32 @@ in the inputs file. Reactions are enabled by setting::
the efficiency).

It is possible to set the maximum and minimum temperature and density for allowing
reactions to occur in a zone using the parameters ``castro.react_T_min``,
``castro.react_T_max``, ``castro.react_rho_min`` and ``castro.react_rho_max``.
reactions to occur in a zone using the parameters:

* ``castro.react_T_min`` and ``castro.react_T_max`` for temperature

* ``castro.react_rho_min`` and ``castro.react_rho_max`` for density

.. index:: castro.disable_shock_burning, USE_SHOCK_VAR

Burning can also be disabled inside shocks. This requires that the code be
compiled with::

USE_SHOCK_VAR = TRUE

in the ``GNUmakefile``. This will allocate storage for a shock flag in the conserved
state array. This flag is computed via a multidimensional shock detection algorithm
that looks for compression (:math:`\nabla \cdot \ub < 0`) along with a pressure jump
in the direction of compression. The runtime parameter::

castro.disable_shock_burning = 1

will skip reactions in a zone where we've detected a shock.

.. note::

Both the compilation with ``USE_SHOCK_VAR = TRUE`` and the runtime parameter
``castro.disable_shock_burning = 1`` are needed to turn off burning in shocks.

Reactions Flowchart
===================
Expand Down
18 changes: 9 additions & 9 deletions Exec/reacting_tests/reacting_convergence/ci-benchmarks/react_converge_128_true_sdc.out
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,16 +21,16 @@
rho_e_sdc_source -6.1984065959e+22 6.3082607681e+21
Temp_sdc_source 0 0
rho_He4_sdc_source -543231.29383 78510.973919
rho_C12_sdc_source -1694.4890575 10.073635142
rho_C12_sdc_source -1694.4890575 10.073635144
rho_O16_sdc_source -0.0098093776861 7.1506350196e-06
rho_Fe56_sdc_source -5.4492579275e-05 7.8513247889e-06
pressure 1.4155320805e+22 4.2608130858e+22
kineng 5.228354476e-13 1.6647276226e+18
kineng 3.6057617076e-14 1.6647276226e+18
soundspeed 212069503.63 257404342.21
Gamma_1 1.5601126452 1.5885713572
MachNumber 6.8192135042e-18 0.0086982771596
MachNumber 1.7908132003e-18 0.0086982771596
entropy 348439780.9 349209883.57
magvort 6.3108872418e-30 0.00018541051835
magvort 0 0.00018541051835
divu -0.1163387912 0.55816306007
eint_E 4.5262357143e+16 6.9937847678e+16
eint_e 4.5262357143e+16 6.9937843728e+16
Expand All @@ -49,11 +49,11 @@
z_velocity 0 0
t_sound_t_enuc 0.00023710316663 0.0195732693
enuc 2.9131490847e+15 4.5102586513e+17
magvel 1.446147223e-09 2067446.1363
radvel -0.00067837194793 2067446.1363
magvel 3.7977686648e-10 2067446.1363
radvel -0.00067839201193 2067446.1363
circvel 0 11.820144682
magmom 0.00072307361144 1.6422547233e+12
angular_momentum_x -0 -0
magmom 0.00018988843323 1.6422547233e+12
angular_momentum_x 0 0
angular_momentum_y 0 0
angular_momentum_z -1.2410862734e+14 1.2410862747e+14
angular_momentum_z -1.241086281e+14 1.2410862748e+14

18 changes: 0 additions & 18 deletions Source/radiation/CastroRad_1d.f90
View file

This file was deleted.

30 changes: 0 additions & 30 deletions Source/radiation/CastroRad_2d.f90
View file

This file was deleted.

18 changes: 0 additions & 18 deletions Source/radiation/CastroRad_3d.f90
View file

This file was deleted.

216 changes: 216 additions & 0 deletions Source/radiation/HABEC.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -367,6 +367,222 @@ namespace HABEC
}
});
}

AMREX_INLINE
void hdterm (Array4<Real> const dterm,
Array4<Real> const er,
const Box& reg,
int cdir, int bct, Real bcl,
Array4<Real const> const bcval,
Array4<int const> const mask,
Array4<Real const> const d,
const Real* dx)
{
Real h;

// Index shift for whichever edge we're checking.
// Negative means we're looking at the lo edge,
// positive means we're looking at the hi edge.
// The first group is for cell centers, the second
// group is for cell edges.

int icp = 0;
int jcp = 0;
int kcp = 0;

int iep = 0;
int jep = 0;
int kep = 0;

#if AMREX_SPACEDIM == 1
if (cdir == 0) {
h = dx[0];
icp = -1;
}
else if (cdir == 1) {
h = dx[0];
icp = 1;
iep = 1;
}
#elif AMREX_SPACEDIM == 2
if (cdir == 0) {
h = dx[0];
icp = -1;
}
else if (cdir == 2) {
h = dx[0];
icp = 1;
iep = 1;
}
else if (cdir == 1) {
h = dx[1];
jcp = -1;
}
else if (cdir == 3) {
h = dx[1];
jcp = 1;
jep = 1;
}
#else
if (cdir == 0) {
h = dx[0];
icp = -1;
}
else if (cdir == 3) {
h = dx[0];
icp = 1;
iep = 1;
}
else if (cdir == 1) {
h = dx[1];
jcp = -1;
}
else if (cdir == 4) {
h = dx[1];
jcp = 1;
jep = 1;
}
else if (cdir == 2) {
h = dx[2];
kcp = -1;
}
else if (cdir == 5) {
h = dx[2];
kcp = 1;
kep = 1;
}
#endif

amrex::LoopOnCpu(reg, [=] (int i, int j, int k) noexcept
{
if (mask.contains(i+icp,j+jcp,k+kcp)) {
if (mask(i+icp,j+jcp,k+kcp) > 0) {
if (bct == LO_DIRICHLET) {
Real d_sign = 1.0_rt;
if (iep != 0 || jep != 0 || kep != 0) {
// right edge
d_sign = -1.0_rt;
}
dterm(i+iep,j+jep,k+kep) = d(i+iep,j+jep,k+kep) * d_sign * (er(i,j,k) - bcval(i+icp,j+jcp,k+kcp)) / (0.5_rt * h + bcl);
}
}
}
});
}

AMREX_INLINE
void hdterm3 (Array4<Real> const dterm,
Array4<Real const> const er,
const Box& reg,
int cdir, int bctype,
Array4<int const> const tf,
Real bcl,
Array4<Real const> const bcval,
Array4<int const> const mask,
Array4<Real const> const d,
const Real* const dx)
{
Real h;

// Index shift for whichever edge we're checking.
// Negative means we're looking at the lo edge,
// positive means we're looking at the hi edge.
// The first group is for cell centers, the second
// group is for cell edges.

int icp = 0;
int jcp = 0;
int kcp = 0;

int iep = 0;
int jep = 0;
int kep = 0;

#if AMREX_SPACEDIM == 1
if (cdir == 0) {
h = dx[0];
icp = -1;
}
else if (cdir == 1) {
h = dx[0];
icp = 1;
iep = 1;
}
#elif AMREX_SPACEDIM == 2
if (cdir == 0) {
h = dx[0];
icp = -1;
}
else if (cdir == 2) {
h = dx[0];
icp = 1;
iep = 1;
}
else if (cdir == 1) {
h = dx[1];
jcp = -1;
}
else if (cdir == 3) {
h = dx[1];
jcp = 1;
jep = 1;
}
#else
if (cdir == 0) {
h = dx[0];
icp = -1;
}
else if (cdir == 3) {
h = dx[0];
icp = 1;
iep = 1;
}
else if (cdir == 1) {
h = dx[1];
jcp = -1;
}
else if (cdir == 4) {
h = dx[1];
jcp = 1;
jep = 1;
}
else if (cdir == 2) {
h = dx[2];
kcp = -1;
}
else if (cdir == 5) {
h = dx[2];
kcp = 1;
kep = 1;
}
#endif

amrex::LoopOnCpu(reg, [=] (int i, int j, int k) noexcept
{
if (mask.contains(i+icp,j+jcp,k+kcp)) {
if (mask(i+icp,j+jcp,k+kcp) > 0) {
int bct;
if (bctype == -1) {
bct = tf(i+icp,j+jcp,k+kcp);
}
else {
bct = bctype;
}
if (bct == LO_DIRICHLET) {
Real d_sign = 1.0_rt;
if (iep != 0 || jep != 0 || kep != 0) {
// right edge
d_sign = -1.0_rt;
}
dterm(i+iep,j+jep,k+kep) = d(i+iep,j+jep,k+kep) * d_sign * (er(i,j,k) - bcval(i+icp,j+jcp,k+kcp)) / (0.5_rt * h + bcl);
}
else if (bct == LO_NEUMANN && bcval(i+icp,j+jcp,k+kcp) == 0.0_rt) {
dterm(i+iep,j+jep,k+kep) = 0.0_rt;
}
}
}
});
}
} // namespace HABEC

#endif
Loading

0 comments on commit 9328a6d

Please sign in to comment.