Skip to content

Commit

Permalink
don't allocate memory in prim state q for passives (#2179)
Browse files Browse the repository at this point in the history
We can create the primitive variable passives as needed in the tracing routines
from the conserved state, saying on memory (with ghost cells) in the primitive
variable array q
  • Loading branch information
zingale authored Jul 1, 2022
1 parent da0f231 commit d68821a
Show file tree
Hide file tree
Showing 14 changed files with 183 additions and 74 deletions.
30 changes: 15 additions & 15 deletions Source/driver/_variables
Original file line number Diff line number Diff line change
Expand Up @@ -24,24 +24,24 @@

# the primitive variable state
@set: primitive NQ
density QRHO NQSRC 1 None
x-velocity QU NQSRC 1 None
y-velocity QV NQSRC 1 None
z-velocity QW NQSRC 1 None
gamma_c QGC NQSRC 1 TRUE_SDC
pressure QPRES NQSRC 1 None
rho-e QREINT NQSRC 1 None
B-x QMAGX NQSRC 1 MHD
B-y QMAGY NQSRC 1 MHD
B-z QMAGZ NQSRC 1 MHD
temperature QTEMP NQSRC 1 None
density QRHO [NQTHERM, NQSRC] 1 None
x-velocity QU [NQTHERM, NQSRC] 1 None
y-velocity QV [NQTHERM, NQSRC] 1 None
z-velocity QW [NQTHERM, NQSRC] 1 None
gamma_c QGC [NQTHERM, NQSRC] 1 TRUE_SDC
pressure QPRES [NQTHERM, NQSRC] 1 None
rho-e QREINT [NQTHERM, NQSRC] 1 None
B-x QMAGX [NQTHERM, NQSRC] 1 MHD
B-y QMAGY [NQTHERM, NQSRC] 1 MHD
B-z QMAGZ [NQTHERM, NQSRC] 1 MHD
temperature QTEMP [NQTHERM, NQSRC] 1 None
total-pressure QPTOT NQTHERM 1 RADIATION
total-pressure QPTOT NQTHERM 1 MHD
total-reint QREITOT NQTHERM 1 RADIATION
radiation QRAD NQTHERM NGROUPS RADIATION
advected QFA None NumAdv None
species QFS None NumSpec None
auxiliary QFX None NumAux None
total-pressure QPTOT None 1 RADIATION
total-pressure QPTOT None 1 MHD
total-reint QREITOT None 1 RADIATION
radiation QRAD None NGROUPS RADIATION


# the auxiliary quantities
Expand Down
52 changes: 32 additions & 20 deletions Source/driver/set_variables.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,18 +7,16 @@
The indices are all 0-based.
For the adds-to column, we can take a list of counters or tuples, e.g.,
[N1, (N2, DEFINE)]
where we only add to N2 is DEFINE is defined as a preprocessor variable.
"""

import argparse
import os
import re


def split_pair(pair_string):
"""given an option of the form "(val1, val2)", split it into val1 and
val2"""
return pair_string.replace("(", "").replace(")", "").replace(" ", "").split(",")

class Index:
"""an index that we want to set"""

Expand Down Expand Up @@ -125,35 +123,44 @@ def doit(variables_file, odir, defines, nadv):
else:

# this splits the line into separate fields. A field is a
# single word or a pair in parentheses like "(a, b)"
fields = re.findall(r'[\w\"\+\.-]+|\([\w+\.-]+\s*,\s*[\w\+\.-]+\)', line)
# single word or, for the "also adds to" section, a group in []
fields = re.findall(r'[\w-]+|\[.*?\]', line)

name = fields[0]
var = fields[1]
adds_to = fields[2]
adds_to_temp = fields[2]
count = fields[3]
ifdef = fields[4]

# we may be fed a pair of the form (SET, DEFINE),
# in which case we only add to SET if we define
# DEFINE
if adds_to.startswith("("):
add_set, define = split_pair(adds_to)
if not define in defines:
if adds_to_temp.startswith("["):
# now split it into entities that are either single
# (counter, ifdef) or just (counter)
subfields = re.findall(r'[\w\"\+\.-]+|\([\w+\.-]+\s*,\s*[\w\+\.-]+\)', adds_to_temp)

# loop over the subfields and build a list of tuples
adds_to = []
for s in subfields:
if s.startswith("("):
counter, define = re.findall(r'[\w]+', s)
if define in defines:
adds_to.append(counter)
else:
adds_to.append(s)
else:
if adds_to_temp == "None":
adds_to = None
else:
adds_to = add_set

if adds_to == "None":
adds_to = None
adds_to = [adds_to_temp]

# only recognize the index if we defined any required preprocessor variable
if ifdef == "None" or ifdef in defines:
indices.append(Index(name, var, default_group=default_group,
iset=current_set, also_adds_to=adds_to,
count=count))


# find the set of set names
unique_sets = {q.iset for q in indices}

Expand All @@ -170,7 +177,13 @@ def doit(variables_file, odir, defines, nadv):
set_indices = [q for q in indices if q.iset == s]

# these indices may also add to other counters
adds_to = {q.adds_to for q in set_indices if q.adds_to is not None}
adds_to = []
for q in set_indices:
if q.adds_to is not None:
for v in q.adds_to:
adds_to.append(v)

adds_to = set(adds_to)

# initialize the counters
counter_main = Counter(default_set[s])
Expand All @@ -187,12 +200,11 @@ def doit(variables_file, odir, defines, nadv):

# increment the counters
counter_main.add_index(i)
if i.adds_to is not None:
if i.adds_to:
for ca in counter_adds:
if ca.name == i.adds_to:
if ca.name in i.adds_to:
ca.add_index(i)


# store the counters for later writing
all_counters += [counter_main]
all_counters += counter_adds
Expand Down
24 changes: 15 additions & 9 deletions Source/hydro/Castro_ctu.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,8 @@ Castro::consup_hydro(const Box& bx,

void
Castro::ctu_ppm_states(const Box& bx, const Box& vbx,
Array4<Real const> const& U_arr,
Array4<Real const> const& rho_inv_arr,
Array4<Real const> const& q_arr,
Array4<Real const> const& qaux_arr,
Array4<Real const> const& srcQ,
Expand Down Expand Up @@ -116,7 +118,7 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx,
if (idir == 0) {
trace_ppm(bx,
idir,
q_arr, qaux_arr, srcQ,
U_arr, rho_inv_arr, q_arr, qaux_arr, srcQ,
qxm, qxp,
#if AMREX_SPACEDIM <= 2
dloga,
Expand All @@ -129,7 +131,7 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx,
} else if (idir == 1) {
trace_ppm(bx,
idir,
q_arr, qaux_arr, srcQ,
U_arr, rho_inv_arr, q_arr, qaux_arr, srcQ,
qym, qyp,
#if AMREX_SPACEDIM <= 2
dloga,
Expand All @@ -144,7 +146,7 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx,
} else {
trace_ppm(bx,
idir,
q_arr, qaux_arr, srcQ,
U_arr, rho_inv_arr, q_arr, qaux_arr, srcQ,
qzm, qzp,
vbx, dt);

Expand All @@ -158,6 +160,8 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx,
#ifdef RADIATION
void
Castro::ctu_ppm_rad_states(const Box& bx, const Box& vbx,
Array4<Real const> const& U_arr,
Array4<Real const> const& rho_inv_arr,
Array4<Real const> const& q_arr,
Array4<Real const> const& qaux_arr,
Array4<Real const> const& srcQ,
Expand All @@ -183,7 +187,7 @@ Castro::ctu_ppm_rad_states(const Box& bx, const Box& vbx,

trace_ppm_rad(bx,
idir,
q_arr, qaux_arr, srcQ,
U_arr, rho_inv_arr, q_arr, qaux_arr, srcQ,
qxm, qxp,
#if AMREX_SPACEDIM <= 2
dloga,
Expand All @@ -196,7 +200,7 @@ Castro::ctu_ppm_rad_states(const Box& bx, const Box& vbx,
} else if (idir == 1) {
trace_ppm_rad(bx,
idir,
q_arr, qaux_arr, srcQ,
U_arr, rho_inv_arr, q_arr, qaux_arr, srcQ,
qym, qyp,
#if AMREX_SPACEDIM <= 2
dloga,
Expand All @@ -211,7 +215,7 @@ Castro::ctu_ppm_rad_states(const Box& bx, const Box& vbx,
} else {
trace_ppm_rad(bx,
idir,
q_arr, qaux_arr, srcQ,
U_arr, rho_inv_arr, q_arr, qaux_arr, srcQ,
qzm, qzp,
vbx, dt);

Expand All @@ -226,6 +230,8 @@ Castro::ctu_ppm_rad_states(const Box& bx, const Box& vbx,

void
Castro::ctu_plm_states(const Box& bx, const Box& vbx,
Array4<Real const> const& U_arr,
Array4<Real const> const& rho_inv_arr,
Array4<Real const> const& q_arr,
Array4<Real const> const& qaux_arr,
Array4<Real const> const& srcQ,
Expand Down Expand Up @@ -267,7 +273,7 @@ Castro::ctu_plm_states(const Box& bx, const Box& vbx,

if (idir == 0) {
trace_plm(bx, 0,
q_arr, qaux_arr,
U_arr, rho_inv_arr, q_arr, qaux_arr,
qxm, qxp,
#if AMREX_SPACEDIM < 3
dloga,
Expand All @@ -279,7 +285,7 @@ Castro::ctu_plm_states(const Box& bx, const Box& vbx,
#if AMREX_SPACEDIM >= 2
} else if (idir == 1) {
trace_plm(bx, 1,
q_arr, qaux_arr,
U_arr, rho_inv_arr, q_arr, qaux_arr,
qym, qyp,
#if AMREX_SPACEDIM < 3
dloga,
Expand All @@ -292,7 +298,7 @@ Castro::ctu_plm_states(const Box& bx, const Box& vbx,
#if AMREX_SPACEDIM == 3
} else {
trace_plm(bx, 2,
q_arr, qaux_arr,
U_arr, rho_inv_arr, q_arr, qaux_arr,
qzm, qzp,
srcQ, vbx, dt);

Expand Down
30 changes: 26 additions & 4 deletions Source/hydro/Castro_ctu_hydro.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -134,6 +134,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)

FArrayBox shk;
FArrayBox q, qaux;
FArrayBox rho_inv;
FArrayBox src_q;
FArrayBox qxm, qxp;
#if AMREX_SPACEDIM >= 2
Expand Down Expand Up @@ -182,7 +183,13 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)
const Box& qbx = amrex::grow(bx, NUM_GROW);
const Box& qbx3 = amrex::grow(bx, 3);

#ifdef RADIATION
q.resize(qbx, NQ);
#else
// note: we won't store the passives in q, so we'll compute their
// primitive versions on demand as needed
q.resize(qbx, NQTHERM);
#endif
Elixir elix_q = q.elixir();
fab_size += q.nBytes();
Array4<Real> const q_arr = q.array();
Expand All @@ -192,7 +199,20 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)
fab_size += qaux.nBytes();
Array4<Real> const qaux_arr = qaux.array();

ctoprim(qbx, time, Sborder.array(mfi),
Array4<Real const> const U_old_arr = Sborder.array(mfi);

rho_inv.resize(qbx3, 1);
Elixir elix_rho_inv = rho_inv.elixir();
fab_size += rho_inv.nBytes();
Array4<Real> const rho_inv_arr = rho_inv.array();

amrex::ParallelFor(qbx3,
[=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
{
rho_inv_arr(i,j,k) = 1.0 / U_old_arr(i,j,k,URHO);
});

ctoprim(qbx, time, U_old_arr,
#ifdef RADIATION
Erborder.array(mfi), lamborder.array(mfi),
#endif
Expand Down Expand Up @@ -264,7 +284,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)
amrex::ParallelFor(qbx3,
[=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
{
hydro::src_to_prim(i, j, k, dt, q_arr, old_src_arr, src_corr_arr, src_q_arr);
hydro::src_to_prim(i, j, k, dt, U_old_arr, q_arr, old_src_arr, src_corr_arr, src_q_arr);
});

// work on the interface states
Expand Down Expand Up @@ -311,6 +331,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)
if (ppm_type == 0) {

ctu_plm_states(obx, bx,
U_old_arr, rho_inv_arr,
q_arr,
qaux_arr,
src_q_arr,
Expand All @@ -330,6 +351,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)

#ifdef RADIATION
ctu_ppm_rad_states(obx, bx,
U_old_arr, rho_inv_arr,
q_arr, qaux_arr, src_q_arr,
qxm_arr, qxp_arr,
#if AMREX_SPACEDIM >= 2
Expand All @@ -345,6 +367,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)
#else

ctu_ppm_states(obx, bx,
U_old_arr, rho_inv_arr,
q_arr, qaux_arr, src_q_arr,
qxm_arr, qxp_arr,
#if AMREX_SPACEDIM >= 2
Expand Down Expand Up @@ -1141,7 +1164,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)
const Box& nbx = amrex::surroundingNodes(bx, idir);

Array4<Real> const flux_arr = (flux[idir]).array();
Array4<Real const> const uin_arr = Sborder.array(mfi);

// Zero out shock and temp fluxes -- these are physically meaningless here
amrex::ParallelFor(nbx,
Expand All @@ -1153,7 +1175,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)
#endif
});

apply_av(nbx, idir, div_arr, uin_arr, flux_arr);
apply_av(nbx, idir, div_arr, U_old_arr, flux_arr);

#ifdef RADIATION
Array4<Real> const rad_flux_arr = (rad_flux[idir]).array();
Expand Down
Loading

0 comments on commit d68821a

Please sign in to comment.