diff --git a/.github/workflows/check-ifdefs.yml b/.github/workflows/check-ifdefs.yml new file mode 100644 index 0000000000..5558864238 --- /dev/null +++ b/.github/workflows/check-ifdefs.yml @@ -0,0 +1,37 @@ +name: check ifdefs + +on: + push: + branches: + - development + - main + pull_request: + branches: + - development + +jobs: + check-ifdefs: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v3 + with: + fetch-depth: 0 + + - name: Setup Python + uses: actions/setup-python@v4 + with: + python-version: '3.10' + + - name: Cache pip + uses: actions/cache@v3 + with: + # this path is specific to Ubuntu + path: ~/.cache/pip + key: ${{ runner.os }}-pip-${{ hashFiles('**/requirements.txt') }} + restore-keys: | + ${{ runner.os }}-pip- + + - name: Run check-ifdefs + run: | + python .github/workflows/check_ifdefs.py .github/workflows/good_defines.txt diff --git a/.github/workflows/check_ifdefs.py b/.github/workflows/check_ifdefs.py new file mode 100644 index 0000000000..2c66ebc96e --- /dev/null +++ b/.github/workflows/check_ifdefs.py @@ -0,0 +1,109 @@ +#!/bin/env python + +import re +import sys + +from pathlib import Path + + +def find_source_files(): + p = Path("./") + files = list(p.glob(r"**/*.cpp")) + files += list(p.glob(r"**/*.H")) + return files + +def check_file(filename): + + # this is a general check to see if we should further examine a line + if_general_re = re.compile(r"^(?:#if|#elif)", re.IGNORECASE|re.DOTALL) + + # this checks something of the form "#ifdef NAME" + ifdef_re = re.compile(r"^#if[n]*def\s+([a-z_0-9]+)", re.IGNORECASE|re.DOTALL) + + # this checks something of the form + # #if (NAME == X) + if_re = re.compile(r"^(?:#if|#elif)\s+[\(]?([a-z_0-9]+)", re.IGNORECASE|re.DOTALL) + + # together these check something of the form + # #if defined(NAME1) || !defined(NAME2) + if_defined_re1 = re.compile(r"^(?:#if|#elif)\s+[!]?(?:defined)", re.IGNORECASE|re.DOTALL) + if_defined_re2 = re.compile(r"[!]?(?:defined)\s*[\(]?([a-z_0-9]+)[\)]?", re.IGNORECASE|re.DOTALL) + + ierr = 0 + defines = [] + + with open(filename) as cf: + for line in cf: + if if_general_re.search(line): + + # check each of the patterns + if if_defined_re1.search(line): + g = if_defined_re2.findall(line) + defines += g + continue + + if g := ifdef_re.search(line): + defines.append(g.group(1)) + continue + + if g := if_re.search(line): + defines.append(g.group(1)) + continue + + + # if we made it here, then we didn't handle things + ierr = 1 + print(f"unhandled, file: {filename} | {line}") + + return ierr, set(defines) + +if __name__ == "__main__": + + good_defines_file = sys.argv[1] + + # read in the list of good defines + + good_defines = [] + with open(good_defines_file) as gd: + for line in gd: + good_defines.append(line.strip()) + + all_defines = [] + total_errors = 0 + for f in find_source_files(): + if "tmp_build_dir" in f.parts: + # skip generated files + continue + ierr, defines = check_file(f) + all_defines += defines + total_errors += ierr + + # remove any header guards + + defines = [] + for d in all_defines: + if d.endswith("_H") or d.endswith("_H_"): + continue + if len(d) == 1 and d.isdigit(): + continue + defines.append(d) + + defines = sorted(set(defines)) + + print("found defines:") + for d in defines: + print(d) + + # now check to make sure that all the defines we found are okay + + invalid = [] + for d in defines: + if d not in good_defines: + invalid.append(d) + + if invalid or total_errors > 0: + if invalid: + print("\ninvalid defines:") + for bad in invalid: + print(bad) + sys.exit(1) diff --git a/.github/workflows/good_defines.txt b/.github/workflows/good_defines.txt new file mode 100644 index 0000000000..5c72004163 --- /dev/null +++ b/.github/workflows/good_defines.txt @@ -0,0 +1,46 @@ +AMREX_DEBUG +AMREX_PARTICLES +AMREX_SPACEDIM +AMREX_USE_CUDA +AMREX_USE_GPU +AMREX_USE_OMP +AUX_THERMO +AUX_UPDATE +BL_FORT_USE_LOWERCASE +BL_FORT_USE_UNDERSCORE +BL_FORT_USE_UPPERCASE +BL_LANG_FORT +BL_LAZY +BL_USE_SETBUF +CXX_MODEL_PARSER +DIFFUSION +DO_PROBLEM_POST_INIT +DO_PROBLEM_POST_RESTART +DO_PROBLEM_POST_SIMULATION +DO_PROBLEM_POST_TIMESTEP +GRAVITY +GR_GRAV +HYBRID_MOMENTUM +INTEGRATOR +MAESTRO_INIT +MHD +NAUX_NET +NEW_NETWORK_IMPLEMENTATION +NGROUPS +NONAKA_PLOT +NSE +NSE_NET +NSE_TABLE +RADIATION +RAD_INTERP +REACTIONS +ROTATION +SCREENING +SHOCK_VAR +SIMPLIFIED_SDC +SPONGE +STRANG +TRUE_SDC +WIN32 +_OPENMP +__cplusplus diff --git a/Diagnostics/Sedov/GNUmakefile b/Diagnostics/Sedov/GNUmakefile index ccc8ba90a2..d9b9c53ec4 100644 --- a/Diagnostics/Sedov/GNUmakefile +++ b/Diagnostics/Sedov/GNUmakefile @@ -25,9 +25,9 @@ Blocs := . CASTRO_HOME := ../.. -INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Extern/amrdata -include $(AMREX_HOME)/Src/Extern/amrdata/Make.package -vpathdir += $(AMREX_HOME)/Src/Extern/amrdata +#INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Extern/amrdata +#include $(AMREX_HOME)/Src/Extern/amrdata/Make.package +#vpathdir += $(AMREX_HOME)/Src/Extern/amrdata include $(CASTRO_HOME)/Exec/Make.Castro diff --git a/Diagnostics/Sedov/Make.package b/Diagnostics/Sedov/Make.package index 98a9ef378e..e69de29bb2 100644 --- a/Diagnostics/Sedov/Make.package +++ b/Diagnostics/Sedov/Make.package @@ -1,2 +0,0 @@ -ca_f90EXE_sources += sedov_util.f90 -CEXE_headers += Sedov_F.H diff --git a/Diagnostics/Sedov/Sedov_F.H b/Diagnostics/Sedov/Sedov_F.H deleted file mode 100644 index 75d1a439a8..0000000000 --- a/Diagnostics/Sedov/Sedov_F.H +++ /dev/null @@ -1,78 +0,0 @@ -#ifndef _Sedov_F_H_ -#define _Sedov_F_H_ -#include - - -#ifdef __cplusplus -extern "C" -{ -#endif - -void fextract1d(const int* lo, const int* hi, - BL_FORT_FAB_ARG_3D(p),const int* nc_p, - const int nbins, - amrex::Real* dens_bin, amrex::Real* vel_bin, - amrex::Real* pres_bin, amrex::Real* e_bin, - int* imask, const int mask_size, const int r1, - const int dens_comp, const int xmom_comp, - const int pres_comp, const int rhoe_comp); - -void fextract2d_cyl(const int* lo, const int* hi, - BL_FORT_FAB_ARG_3D(p),const int* nc_p, - const int nbins, - amrex::Real* dens_bin, amrex::Real* vel_bin, - amrex::Real* pres_bin, amrex::Real* e_bin, - int* ncount, - int* imask, const int mask_size, const int r1, - const int dens_comp, const int xmom_comp, - const int ymom_comp, - const int pres_comp, const int rhoe_comp, - const amrex::Real dx_fine, const amrex::Real* dx, - const amrex::Real xctr, const amrex::Real yctr); - -void fextract2d_sph(const int* lo, const int* hi, - BL_FORT_FAB_ARG_3D(p),const int* nc_p, - const int nbins, - amrex::Real* dens_bin, amrex::Real* vel_bin, - amrex::Real* pres_bin, amrex::Real* e_bin, - amrex::Real* volcount, - int* imask, const int mask_size, const int r1, - const int dens_comp, const int xmom_comp, - const int ymom_comp, - const int pres_comp, const int rhoe_comp, - const amrex::Real dx_fine, const amrex::Real* dx, - const amrex::Real xctr, const amrex::Real yctr); - -void fextract3d_cyl(const int* lo, const int* hi, - BL_FORT_FAB_ARG_3D(p),const int* nc_p, - const int nbins, - amrex::Real* dens_bin, amrex::Real* vel_bin, - amrex::Real* pres_bin, - int* ncount, - int* imask, const int mask_size, const int r1, - const int dens_comp, const int xmom_comp, - const int ymom_comp, const int zmom_comp, - const int pres_comp, - const amrex::Real dx_fine, const amrex::Real* dx, - const amrex::Real xctr, const amrex::Real yctr); - -void fextract3d_sph(const int* lo, const int* hi, - BL_FORT_FAB_ARG_3D(p),const int* nc_p, - const int nbins, - amrex::Real* dens_bin, amrex::Real* vel_bin, - amrex::Real* pres_bin, amrex::Real* e_bin, - int* ncount, - - int* imask, const int mask_size, const int r1, - const int dens_comp, const int xmom_comp, - const int ymom_comp, const int zmom_comp, - const int pres_comp, const int rhoe_comp, - const amrex::Real dx_fine, const amrex::Real* dx, - const amrex::Real xctr, - const amrex::Real yctr, const amrex::Real zctr); - -#ifdef __cplusplus -} -#endif - -#endif diff --git a/Diagnostics/Sedov/main.cpp b/Diagnostics/Sedov/main.cpp index 10d332634e..e40751b6c9 100644 --- a/Diagnostics/Sedov/main.cpp +++ b/Diagnostics/Sedov/main.cpp @@ -5,8 +5,10 @@ #include // #include #include -#include -#include +#include +#include +#include +#include using namespace amrex; @@ -16,360 +18,443 @@ std::string inputs_name = ""; // Prototypes // void GetInputArgs (const int argc, char** argv, - string& pltfile, string& slcfile, + std::string& pltfile, std::string& slcfile, bool& sphr); -string GetVarFromJobInfo (const string pltfile, const string varname); +std::string GetVarFromJobInfo (const std::string pltfile, const std::string varname); -Vector GetCenter (const string pltfile); +Vector GetCenter (const std::string pltfile); void PrintHelp (); -int main(int argc, char* argv[]) -{ +std::pair get_coord_info(const Array& p, + const Vector& center, + const Array& dx_level, + const int coord, const bool sphr) { - amrex::Initialize(argc, argv, false); - // timer for profiling - BL_PROFILE_VAR("main()", pmain); + // compute radial coord and index - // Input arguments - string pltfile, slcfile; - bool sphr = false; + Real r_zone{0.0_rt}; + Real vol = std::numeric_limits::quiet_NaN(); - GetInputArgs (argc, argv, pltfile, slcfile, sphr); +#if AMREX_SPACEDIM == 1 + // 1-d spherical geometry / spherical Sedov explosion - auto center = GetCenter(pltfile); - double xctr = center[0]; - double yctr = center[1]; - double zctr = center[2]; - // Start dataservices - DataServices::SetBatchMode(); + AMREX_ASSERT(coord == 2); - // Define the type of file - Amrvis::FileType fileType(Amrvis::NEWPLT); - DataServices dataServices (pltfile, fileType); + r_zone = p[0] - center[0]; + Real r_r = problo[0]+static_cast(i+1)*dx_level[0]; + Real r_l = problo[0]+static_cast(i)*dx_level[0]; + vol = (4.0_rt/3.0_rt) * M_PI * dx_level[0] * + (r_r*r_r + r_l*r_r + r_l*r_l); - if (!dataServices.AmrDataOk()) - DataServices::Dispatch(DataServices::ExitRequest, NULL); +#elif AMREX_SPACEDIM == 2 + if (sphr) { + // 2-d axisymmetric geometry / spherical Sedov explosion - // get data from plot file - AmrData& data = dataServices.AmrDataRef(); + AMREX_ASSERT(coord == 1); - int finestLevel = data.FinestLevel(); + // axisymmetric V = pi (r_r**2 - r_l**2) * dz + // = pi dr * dz * (r_r + r_l) + // = 2 pi r dr dz - // get variable names - const Vector& varNames = data.PlotVarNames(); + r_zone = std::sqrt((p[0] - center[0]) * (p[0] - center[0]) + + (p[1] - center[1]) * (p[1] - center[1])); + vol = 2 * M_PI * p[0] * dx_level[0] * dx_level[1]; - // get the index bounds and dx. - Box domain = data.ProbDomain()[finestLevel]; - Vector dx = data.CellSize(finestLevel); - const Vector& problo = data.ProbLo(); - const Vector& probhi = data.ProbHi(); + } else { + // 2-d Cartesian geometry / cylindrical Sedov explosion - Vector rr = data.RefRatio(); + AMREX_ASSERT(coord == 0); + + r_zone = std::sqrt((p[0] - center[0]) * (p[0] - center[0]) + + (p[1] - center[1]) * (p[1] - center[1])); + vol = dx_level[0] * dx_level[1]; + + } - // compute the size of the radially-binned array -- we'll do it to - // the furtherest corner of the domain -#if (AMREX_SPACEDIM == 1) - double maxdist = fabs(probhi[0] - problo[0]); -#elif (AMREX_SPACEDIM == 2) - double x_maxdist = max(fabs(probhi[0] - xctr), fabs(problo[0] - xctr)); - double y_maxdist = max(fabs(probhi[1] - yctr), fabs(problo[1] - yctr)); - double maxdist = sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist); #else - double x_maxdist = max(fabs(probhi[0] - xctr), fabs(problo[0] - xctr)); - double y_maxdist = max(fabs(probhi[1] - yctr), fabs(problo[1] - yctr)); - double z_maxdist = max(fabs(probhi[2] - zctr), fabs(problo[2] - zctr)); - double maxdist = sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist + - z_maxdist*z_maxdist); -#endif - double dx_fine = *(std::min_element(dx.begin(), dx.end())); + AMREX_ASSERT(coord == 0); - int nbins = int(maxdist / dx_fine); + vol = dx_level[0] * dx_level[1] * dx_level[2]; - // radial coordinate - Vector r(nbins); + if (sphr) { + // 3-d Cartesian geometry / spherical Sedov explosion - for (auto i = 0; i < nbins; i++) - r[i] = (i + 0.5) * dx_fine; + r_zone = std::sqrt((p[0] - center[0]) * (p[0] - center[0]) + + (p[1] - center[1]) * (p[1] - center[1]) + + (p[2] - center[2]) * (p[2] - center[2])); - // find variable indices - auto dens_comp = data.StateNumber("density"); - auto xmom_comp = data.StateNumber("xmom"); -#if (AMREX_SPACEDIM >= 2) - auto ymom_comp = data.StateNumber("ymom"); -#endif -#if (AMREX_SPACEDIM == 3) - auto zmom_comp = data.StateNumber("zmom"); -#endif - auto pres_comp = data.StateNumber("pressure"); - auto rhoe_comp = data.StateNumber("rho_e"); + } else { + // 3-d Cartesian geometry / cylindrical Sedov explosion - if (dens_comp < 0 || xmom_comp < 0 || pres_comp < 0 || rhoe_comp < 0) - Abort("ERROR: variable(s) not found"); + r_zone = std::sqrt((p[0] - center[0]) * (p[0] - center[0]) + + (p[1] - center[1]) * (p[1] - center[1])); -#if (AMREX_SPACEDIM == 3) - if (ymom_comp < 0 || zmom_comp < 0) - Abort("ERROR: variable(s) not found"); + } #endif - // allocate storage for data - Vector dens_bin(nbins, 0.); - Vector vel_bin(nbins, 0.); - Vector pres_bin(nbins, 0.); - Vector e_bin(nbins, 0.); - Vector ncount(nbins, 0); - Vector volcount(nbins, 0); - - // r1 is the factor between the current level grid spacing and the - // FINEST level - auto r1 = 1.0; - - // fill a multifab with the data - Vector fill_comps(data.NComp()); - for (auto i = 0; i < data.NComp(); i++) - fill_comps[i] = i; - - // imask will be set to false if we've already output the data. - // Note, imask is defined in terms of the finest level. As we loop - // over levels, we will compare to the finest level index space to - // determine if we've already output here - int mask_size = domain.length().max(); - Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); - - // loop over the data, starting at the finest grid, and if we haven't - // already stored data in that grid location (according to imask), - // store it. - for (int l = finestLevel; l >= 0; l--) { - - Vector level_dx = data.DxLevel()[l]; - - const BoxArray& ba = data.boxArray(l); - const DistributionMapping& dm = data.DistributionMap(l); - - MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); - data.FillVar(lev_data_mf, l, varNames, fill_comps); - -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); + return {r_zone, vol}; + +} + +int main(int argc, char* argv[]) +{ + + amrex::Initialize(argc, argv, false); + + // timer for profiling + BL_PROFILE_VAR("main()", pmain); + + // Input arguments + std::string pltfile, slcfile; + bool sphr = false; + + GetInputArgs (argc, argv, pltfile, slcfile, sphr); + + auto center = GetCenter(pltfile); + double xctr = center[0]; + double yctr = center[1]; + double zctr = center[2]; + + PlotFileData pf(pltfile); + + int fine_level = pf.finestLevel(); + const int dim = pf.spaceDim(); + + // get the index bounds and dx. + Box domain = pf.probDomain(fine_level); + auto dx = pf.cellSize(fine_level); + auto problo = pf.probLo(); + auto probhi = pf.probHi(); + + // compute the size of the radially-binned array -- we'll do it to + // the furtherest corner of the domain #if (AMREX_SPACEDIM == 1) - fextract1d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, dens_bin.dataPtr(), - vel_bin.dataPtr(), pres_bin.dataPtr(), - e_bin.dataPtr(), imask.dataPtr(), mask_size, r1, - dens_comp, xmom_comp, pres_comp, rhoe_comp); + double maxdist = std::abs(probhi[0] - problo[0]); + #elif (AMREX_SPACEDIM == 2) - if (sphr) { - - fextract2d_sph(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, dens_bin.dataPtr(), - vel_bin.dataPtr(), pres_bin.dataPtr(), - e_bin.dataPtr(), volcount.dataPtr(), - imask.dataPtr(), mask_size, r1, - dens_comp, xmom_comp, ymom_comp, pres_comp, rhoe_comp, - dx_fine, level_dx.dataPtr(), - xctr, yctr); - - } else { - fextract2d_cyl(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, dens_bin.dataPtr(), - vel_bin.dataPtr(), pres_bin.dataPtr(), - e_bin.dataPtr(), ncount.dataPtr(), - imask.dataPtr(), mask_size, r1, - dens_comp, xmom_comp, ymom_comp, pres_comp, rhoe_comp, - dx_fine, level_dx.dataPtr(), - xctr, yctr); - } + double x_maxdist = amrex::max(std::abs(probhi[0] - xctr), + std::abs(problo[0] - xctr)); + double y_maxdist = amrex::max(std::abs(probhi[1] - yctr), + std::abs(problo[1] - yctr)); + double maxdist = std::sqrt(x_maxdist*x_maxdist + + y_maxdist*y_maxdist); + #else - if (sphr) { - fextract3d_sph(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, dens_bin.dataPtr(), - vel_bin.dataPtr(), pres_bin.dataPtr(), - e_bin.dataPtr(), ncount.dataPtr(), - imask.dataPtr(), mask_size, r1, - dens_comp, xmom_comp, ymom_comp, zmom_comp, pres_comp, rhoe_comp, - dx_fine, level_dx.dataPtr(), xctr, yctr, zctr); - - } else { - fextract3d_cyl(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, dens_bin.dataPtr(), - vel_bin.dataPtr(), pres_bin.dataPtr(), - ncount.dataPtr(), - imask.dataPtr(), mask_size, r1, - dens_comp, xmom_comp, ymom_comp, zmom_comp, pres_comp, - dx_fine, level_dx.dataPtr(), xctr, yctr); - } -#endif - } - - // adjust r1 for the next lowest level - if (l != 0) r1 *= rr[l-1]; - } - -#if (AMREX_SPACEDIM >= 2) - - if (AMREX_SPACEDIM == 2 && sphr) { - //normalize - for (int i = 0; i < nbins; i++) { - if (volcount[i] != 0.0) { - dens_bin[i] /= volcount[i]; - vel_bin[i] /= volcount[i]; - pres_bin[i] /= volcount[i]; - e_bin[i] /= volcount[i]; - } - } - } else { - //normalize - for (int i = 0; i < nbins; i++) { - if (ncount[i] != 0) { - dens_bin[i] /= ncount[i]; - vel_bin[i] /= ncount[i]; - pres_bin[i] /= ncount[i]; - e_bin[i] /= ncount[i]; - } - } - } + double x_maxdist = amrex::max(std::abs(probhi[0] - xctr), + std::fabs(problo[0] - xctr)); + double y_maxdist = amrex::max(std::abs(probhi[1] - yctr), + std::abs(problo[1] - yctr)); + double z_maxdist = amrex::max(std::abs(probhi[2] - zctr), + std::abs(problo[2] - zctr)); + double maxdist = std::sqrt(x_maxdist*x_maxdist + + y_maxdist*y_maxdist + + z_maxdist*z_maxdist); #endif - // now open the slicefile and write out the data - std::ofstream slicefile; - slicefile.open(slcfile); - slicefile.setf(std::ios::scientific); - slicefile.precision(12); - const auto w = 24; + double dx_fine = *(std::min_element(dx.begin(), dx.end())); + + int nbins = int(maxdist / dx_fine); + + // radial coordinate + Vector r(nbins); + + for (auto i = 0; i < nbins; i++) + r[i] = (i + 0.5) * dx_fine; + + // find variable indices + const Vector& var_names_pf = pf.varNames(); + + int dens_comp = std::distance(var_names_pf.cbegin(), + std::find(var_names_pf.cbegin(), var_names_pf.cend(), "density")); + + int xmom_comp = std::distance(var_names_pf.cbegin(), + std::find(var_names_pf.cbegin(), var_names_pf.cend(), "xmom")); + + int ymom_comp = std::distance(var_names_pf.cbegin(), + std::find(var_names_pf.cbegin(), var_names_pf.cend(), "ymom")); + + int zmom_comp = std::distance(var_names_pf.cbegin(), + std::find(var_names_pf.cbegin(), var_names_pf.cend(), "zmom")); + + int pres_comp = std::distance(var_names_pf.cbegin(), + std::find(var_names_pf.cbegin(), var_names_pf.cend(), "pressure")); + + int rhoe_comp = std::distance(var_names_pf.cbegin(), + std::find(var_names_pf.cbegin(), var_names_pf.cend(), "rho_e")); + + int coord = pf.coordSys(); + + // allocate storage for data + Vector dens_bin(nbins, 0.); + Vector vel_bin(nbins, 0.); + Vector pres_bin(nbins, 0.); + Vector e_bin(nbins, 0.); + Vector volcount(nbins, 0); + + // we will use a mask that tells us if a zone on the current level + // is covered by data on a finer level. + + for (int ilev = 0; ilev <= fine_level; ++ilev) { + + Array dx_level = pf.cellSize(ilev); + + if (ilev < fine_level) { + IntVect ratio{pf.refRatio(ilev)}; + for (int idim = dim; idim < AMREX_SPACEDIM; ++idim) { + ratio[idim] = 1; + } + const iMultiFab mask = makeFineMask(pf.boxArray(ilev), pf.DistributionMap(ilev), + pf.boxArray(ilev+1), ratio); + + const MultiFab& lev_data_mf = pf.get(ilev); + + for (MFIter mfi(lev_data_mf); mfi.isValid(); ++mfi) { + const Box& bx = mfi.validbox(); + if (bx.ok()) { + const auto& m = mask.array(mfi); + const auto& fab = lev_data_mf.array(mfi); + const auto lo = amrex::lbound(bx); + const auto hi = amrex::ubound(bx); + + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + if (m(i,j,k) == 0) { // not covered by fine + + Array p + = {AMREX_D_DECL(problo[0]+static_cast(i+0.5)*dx_level[0], + problo[1]+static_cast(j+0.5)*dx_level[1], + problo[2]+static_cast(k+0.5)*dx_level[2])}; + + + const auto& [r_zone, vol] = get_coord_info(p, center, dx_level, coord, sphr); + + int index = static_cast(r_zone / dx_fine); + + // add to the bin, weighting by the size + + dens_bin[index] += fab(i,j,k,dens_comp) * vol; + + vel_bin[index] += + std::sqrt(std::pow(fab(i,j,k,xmom_comp), 2) + + std::pow(fab(i,j,k,ymom_comp), 2) + + std::pow(fab(i,j,k,zmom_comp), 2)) / + fab(i,j,k,dens_comp) * vol; + + pres_bin[index] += fab(i,j,k,pres_comp) * vol; + + e_bin[index] += (fab(i,j,k,rhoe_comp) / fab(i,j,k,dens_comp)) * vol; + + volcount[index] += vol; + + } // mask + + } + } + } + + } // bx.ok() + + } // MFIter + + } else { + // this is the finest level + + const MultiFab& lev_data_mf = pf.get(ilev); + + for (MFIter mfi(lev_data_mf); mfi.isValid(); ++mfi) { + const Box& bx = mfi.validbox(); + if (bx.ok()) { + const auto& fab = lev_data_mf.array(mfi); + const auto lo = amrex::lbound(bx); + const auto hi = amrex::ubound(bx); + + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + + Array p + = {AMREX_D_DECL(problo[0]+static_cast(i+0.5)*dx_level[0], + problo[1]+static_cast(j+0.5)*dx_level[1], + problo[2]+static_cast(k+0.5)*dx_level[2])}; + + const auto& [r_zone, vol] = get_coord_info(p, center, dx_level, coord, sphr); + + int index = static_cast(r_zone / dx_fine); + + // add to the bin, weighting by the size + + dens_bin[index] += fab(i,j,k,dens_comp) * vol; + + vel_bin[index] += + std::sqrt(std::pow(fab(i,j,k,xmom_comp), 2) + + std::pow(fab(i,j,k,ymom_comp), 2) + + std::pow(fab(i,j,k,zmom_comp), 2)) / + fab(i,j,k,dens_comp) * vol; + + pres_bin[index] += fab(i,j,k,pres_comp) * vol; + + e_bin[index] += (fab(i,j,k,rhoe_comp) / fab(i,j,k,dens_comp)) * vol; + + volcount[index] += vol; + + } + } + } + + } // bx.ok() + + } // MFIter + + + } + + } // level loop + + + //normalize + for (int i = 0; i < nbins; i++) { + if (volcount[i] != 0.0) { + dens_bin[i] /= volcount[i]; + vel_bin[i] /= volcount[i]; + pres_bin[i] /= volcount[i]; + e_bin[i] /= volcount[i]; + } + } + + // now open the slicefile and write out the data + std::ofstream slicefile; + slicefile.open(slcfile); + slicefile.setf(std::ios::scientific); + slicefile.precision(12); + const auto w = 24; - // write the header - slicefile << "# " << std::setw(w) << "x" << std::setw(w) << "density" << std::setw(w) << "velocity" << std::setw(w) << "pressure" << std::setw(w) << "int. energy" << std::endl; + // write the header + slicefile << "# " << std::setw(w) << "x" << std::setw(w) << "density" << std::setw(w) << "velocity" << std::setw(w) << "pressure" << std::setw(w) << "int. energy" << std::endl; - // write the data in columns - const auto SMALL = 1.e-20; - for (auto i = 0; i < nbins; i++) { - if (fabs(dens_bin[i]) < SMALL) dens_bin[i] = 0.0; - if (fabs( vel_bin[i]) < SMALL) vel_bin[i] = 0.0; - if (fabs(pres_bin[i]) < SMALL) pres_bin[i] = 0.0; - if (fabs( e_bin[i]) < SMALL) e_bin[i] = 0.0; + // write the data in columns + const auto SMALL = 1.e-20; + for (auto i = 0; i < nbins; i++) { + if (fabs(dens_bin[i]) < SMALL) dens_bin[i] = 0.0; + if (fabs( vel_bin[i]) < SMALL) vel_bin[i] = 0.0; + if (fabs(pres_bin[i]) < SMALL) pres_bin[i] = 0.0; + if (fabs( e_bin[i]) < SMALL) e_bin[i] = 0.0; - slicefile << std::setw(w) << r[i] << std::setw(w) << dens_bin[i] << std::setw(w) << vel_bin[i] << std::setw(w) << pres_bin[i] << std::setw(w) << e_bin[i] << std::endl; - } + slicefile << std::setw(w) << r[i] << std::setw(w) << dens_bin[i] << std::setw(w) << vel_bin[i] << std::setw(w) << pres_bin[i] << std::setw(w) << e_bin[i] << std::endl; + } - slicefile.close(); + slicefile.close(); - // destroy timer for profiling - BL_PROFILE_VAR_STOP(pmain); + // destroy timer for profiling + BL_PROFILE_VAR_STOP(pmain); - amrex::Finalize(); + amrex::Finalize(); } // // Parse command line arguments // void GetInputArgs ( const int argc, char** argv, - string& pltfile, string& slcfile, + std::string& pltfile, std::string& slcfile, bool &sphr) { - int i = 1; // skip program name - - while ( i < argc) - { - - if ( !strcmp(argv[i], "-p") || !strcmp(argv[i],"--pltfile") ) - { - pltfile = argv[++i]; - } - else if ( !strcmp(argv[i], "-s") || !strcmp(argv[i],"--slicefile") ) - { - slcfile = argv[++i]; - } - else if ( !strcmp(argv[i],"--sphr") ) - { - sphr = true; - } - else - { - std::cout << "\n\nOption " << argv[i] << " not recognized" << std::endl; - PrintHelp (); - exit ( EXIT_FAILURE ); - } - - // Go to the next parameter name - ++i; - } - - if (pltfile.empty() && slcfile.empty()) - { - PrintHelp(); - Abort("Missing input file"); - } + int i = 1; // skip program name + + while ( i < argc) + { + + if ( !strcmp(argv[i], "-p") || !strcmp(argv[i],"--pltfile") ) + { + pltfile = argv[++i]; + } + else if ( !strcmp(argv[i], "-s") || !strcmp(argv[i],"--slicefile") ) + { + slcfile = argv[++i]; + } + else if ( !strcmp(argv[i],"--sphr") ) + { + sphr = true; + } + else + { + std::cout << "\n\nOption " << argv[i] << " not recognized" << std::endl; + PrintHelp (); + exit ( EXIT_FAILURE ); + } + + // Go to the next parameter name + ++i; + } + + if (pltfile.empty() && slcfile.empty()) + { + PrintHelp(); + Abort("Missing input file"); + } #if (AMREX_SPACEDIM == 1) - Print() << "Extracting slice from 1d problem" << std::endl; + Print() << "Extracting slice from 1d problem" << std::endl; #elif (AMREX_SPACEDIM >= 2) - if (sphr) { - Print() << "Extracting slice from " << AMREX_SPACEDIM << "d spherical problem" << std::endl; - } else { - Print() << "Extracting slice from " << AMREX_SPACEDIM << "d cylindrical problem" << std::endl; - } + if (sphr) { + Print() << "Extracting slice from " << AMREX_SPACEDIM << "d spherical problem" << std::endl; + } else { + Print() << "Extracting slice from " << AMREX_SPACEDIM << "d cylindrical problem" << std::endl; + } #endif - Print() << "\nplotfile = \"" << pltfile << "\"" << std::endl; - Print() << "slicefile = \"" << slcfile << "\"" << std::endl; - Print() << std::endl; + Print() << "\nplotfile = \"" << pltfile << "\"" << std::endl; + Print() << "slicefile = \"" << slcfile << "\"" << std::endl; + Print() << std::endl; } /// /// Gets the variable ``varname`` from the ``job_info`` file and returns as a /// string /// -string GetVarFromJobInfo (const string pltfile, const string varname) { - string filename = pltfile + "/job_info"; - std::regex re("(?:[ \\t]*)" + varname + "\\s*:\\s*(.*)\\s*\\n"); - - std::smatch m; - - std::ifstream jobfile(filename); - if (jobfile.is_open()) { - std::stringstream buf; - buf << jobfile.rdbuf(); - string file_contents = buf.str(); - - if (std::regex_search(file_contents, m, re)) { - return m[1]; - } else { - Print() << "Unable to find " << varname << " in job_info file!" << std::endl; - } - } else { - Print() << "Could not open job_info file!" << std::endl; - } - - return ""; +std::string GetVarFromJobInfo (const std::string pltfile, const std::string varname) { + std::string filename = pltfile + "/job_info"; + std::regex re("(?:[ \\t]*)" + varname + "\\s*:\\s*(.*)\\s*\\n"); + + std::smatch m; + + std::ifstream jobfile(filename); + if (jobfile.is_open()) { + std::stringstream buf; + buf << jobfile.rdbuf(); + std::string file_contents = buf.str(); + + if (std::regex_search(file_contents, m, re)) { + return m[1]; + } else { + Print() << "Unable to find " << varname << " in job_info file!" << std::endl; + } + } else { + Print() << "Could not open job_info file!" << std::endl; + } + + return ""; } // Get the center from the job info file and return as a Real Vector -Vector GetCenter (const string pltfile) { - auto center_str = GetVarFromJobInfo(pltfile, "center"); +Vector GetCenter (const std::string pltfile) { + auto center_str = GetVarFromJobInfo(pltfile, "center"); - // split string - std::istringstream iss {center_str}; - Vector center; + // split string + std::istringstream iss {center_str}; + Vector center; - std::string s; - while (std::getline(iss, s, ',')) - center.push_back(stod(s)); + std::string s; + while (std::getline(iss, s, ',')) + center.push_back(stod(s)); - return center; + return center; } // @@ -377,12 +462,12 @@ Vector GetCenter (const string pltfile) { // void PrintHelp () { - Print() << "\nusage: executable_name args" - << "\nargs [-p|--pltfile] plotfile : plot file directory (required)" - << "\n [-s|--slicefile] slice file : slice file (required)" + Print() << "\nusage: executable_name args" + << "\nargs [-p|--pltfile] plotfile : plot file directory (required)" + << "\n [-s|--slicefile] slice file : slice file (required)" #if AMREX_SPACEDIM >= 2 - << "\n [--sphr] spherical : spherical problem" + << "\n [--sphr] spherical : spherical problem" #endif - << "\n\n" << std::endl; + << "\n\n" << std::endl; } diff --git a/Diagnostics/Sedov/sedov_util.f90 b/Diagnostics/Sedov/sedov_util.f90 deleted file mode 100644 index 33518fec74..0000000000 --- a/Diagnostics/Sedov/sedov_util.f90 +++ /dev/null @@ -1,369 +0,0 @@ -! Process a sedov problem to produce rho, u, p and e as a -! function of r, for comparison to the analytic solution. - -subroutine fextract1d(lo, hi, p, plo, phi, nc_p, nbins, dens_bin, & - vel_bin, pres_bin, e_bin, imask, mask_size, r1,& - dens_comp, xmom_comp, pres_comp, rhoe_comp) bind(C, name='fextract1d') - - use amrex_fort_module, only : rt => amrex_real - use amrex_constants_module - - implicit none - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: plo(3), phi(3), nc_p - integer, intent(in), value :: nbins - real(rt), intent(in) :: p(plo(1):phi(1),plo(2):phi(2),plo(3):phi(3),0:nc_p-1) - real(rt), intent(inout) :: dens_bin(0:nbins-1) - real(rt), intent(inout) :: vel_bin(0:nbins-1) - real(rt), intent(inout) :: pres_bin(0:nbins-1) - real(rt), intent(inout) :: e_bin(0:nbins-1) - integer, intent(inout) :: imask(0:mask_size-1) - integer, intent(in), value :: mask_size, r1, dens_comp, xmom_comp, pres_comp, rhoe_comp - - integer :: i, j, k, index - - ! loop over all of the zones in the patch. Here, we convert - ! the cell-centered indices at the current level into the - ! corresponding RANGE on the finest level, and test if we've - ! stored data in any of those locations. If we haven't then - ! we store this level's data and mark that range as filled. - j = lo(2) - k = lo(3) - do i = lo(1), hi(1) - - if ( any(imask(i*r1:(i+1)*r1-1) .eq. 1) ) then - - index = i * r1 - - dens_bin(index:index+(r1-1)) = p(i,j,k,dens_comp) - - vel_bin(index:index+(r1-1)) = & - abs(p(i,j,k,xmom_comp)) / p(i,j,k,dens_comp) - - pres_bin(index:index+(r1-1)) = p(i,j,k,pres_comp) - - e_bin(index:index+(r1-1)) = & - p(i,j,k,rhoe_comp) / p(i,j,k,dens_comp) - - imask(i*r1:(i+1)*r1-1) = 0 - - end if - - enddo - -end subroutine fextract1d - -subroutine fextract2d_cyl(lo, hi, p, plo, phi, nc_p, nbins, dens_bin, & - vel_bin, pres_bin, e_bin, ncount, imask, mask_size, r1, & - dens_comp, xmom_comp, ymom_comp, pres_comp, rhoe_comp, dx_fine, dx, & - xctr, yctr) & - bind(C, name='fextract2d_cyl') - - use amrex_fort_module, only : rt => amrex_real - use amrex_constants_module - - implicit none - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: plo(3), phi(3), nc_p - integer, intent(in), value :: nbins - real(rt), intent(in) :: p(plo(1):phi(1),plo(2):phi(2),plo(3):phi(3),0:nc_p-1) - real(rt), intent(inout) :: dens_bin(0:nbins-1) - real(rt), intent(inout) :: vel_bin(0:nbins-1) - real(rt), intent(inout) :: pres_bin(0:nbins-1) - real(rt), intent(inout) :: e_bin(0:nbins-1) - integer, intent(inout) :: ncount(0:nbins-1) - integer, intent(in), value :: mask_size - integer, intent(inout) :: imask(0:mask_size-1,0:mask_size-1) - integer, intent(in), value :: r1, dens_comp, xmom_comp, ymom_comp, pres_comp, rhoe_comp - real(rt), intent(in), value :: dx_fine, xctr, yctr - real(rt), intent(in) :: dx(3) - - integer :: i, j, k, index - real(rt) :: xx, yy, r_zone - - ! loop over all of the zones in the patch. Here, we convert - ! the cell-centered indices at the current level into the - ! corresponding RANGE on the finest level, and test if we've - ! stored data in any of those locations. If we haven't then - ! we store this level's data and mark that range as filled. - k = lo(3) - do j = lo(2), hi(2) - yy = (j + HALF)*dx(2) - - do i = lo(1), hi(1) - xx = (i + HALF)*dx(1) - - if ( any(imask(i*r1:(i+1)*r1-1, & - j*r1:(j+1)*r1-1) .eq. 1) ) then - - r_zone = sqrt((xx-xctr)**2 + (yy-yctr)**2) - - index = r_zone/dx_fine - - dens_bin(index) = dens_bin(index) + & - p(i,j,k,dens_comp)*r1**2 - - vel_bin(index) = vel_bin(index) + & - (sqrt(p(i,j,k,xmom_comp)**2 + & - p(i,j,k,ymom_comp)**2)/ & - p(i,j,k,dens_comp))*r1**2 - - pres_bin(index) = pres_bin(index) + & - p(i,j,k,pres_comp)*r1**2 - - e_bin(index) = e_bin(index) + & - (p(i,j,k,rhoe_comp)/p(i,j,k,dens_comp))*r1**2 - - ncount(index) = ncount(index) + r1**2 - - imask(i*r1:(i+1)*r1-1, & - j*r1:(j+1)*r1-1) = 0 - - end if - - enddo - enddo - -end subroutine fextract2d_cyl - -subroutine fextract2d_sph(lo, hi, p, plo, phi, nc_p, nbins, dens_bin, & - vel_bin, pres_bin, e_bin, volcount, imask, mask_size, r1,& - dens_comp, xmom_comp, ymom_comp, pres_comp, rhoe_comp, dx_fine, dx, xctr, yctr) & - bind(C, name='fextract2d_sph') - - use amrex_fort_module, only : rt => amrex_real - use amrex_constants_module - - implicit none - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: plo(3), phi(3), nc_p - integer, intent(in), value :: nbins - real(rt), intent(in) :: p(plo(1):phi(1),plo(2):phi(2),plo(3):phi(3),0:nc_p-1) - real(rt), intent(inout) :: dens_bin(0:nbins-1) - real(rt), intent(inout) :: vel_bin(0:nbins-1) - real(rt), intent(inout) :: pres_bin(0:nbins-1) - real(rt), intent(inout) :: e_bin(0:nbins-1) - real(rt), intent(inout) :: volcount(0:nbins-1) - integer, intent(in), value :: mask_size - integer, intent(inout) :: imask(0:mask_size-1,0:mask_size-1) - integer, intent(in), value :: r1, dens_comp, xmom_comp, ymom_comp, pres_comp, rhoe_comp - real(rt), intent(in), value :: dx_fine, xctr, yctr - real(rt), intent(in) :: dx(3) - - integer :: i, j, k, index - real(rt) :: xx, xl, xr, yy, yl, yr, r_zone, vol, vel - - ! write(*,*) "imask = ", imask - - ! loop over all of the zones in the patch. Here, we convert - ! the cell-centered indices at the current level into the - ! corresponding RANGE on the finest level, and test if we've - ! stored data in any of those locations. If we haven't then - ! we store this level's data and mark that range as filled. - k = lo(3) - do j = lo(2), hi(2) - yy = (dble(j) + HALF)*dx(2) - yl = (dble(j))*dx(2) - yr = (dble(j) + ONE)*dx(2) - - do i = lo(1), hi(1) - xx = (dble(i) + HALF)*dx(1) - xl = (dble(i))*dx(1) - xr = (dble(i) + ONE)*dx(1) - - if ( any(imask(i*r1:(i+1)*r1-1, & - j*r1:(j+1)*r1-1) .eq. 1) ) then - - r_zone = sqrt((xx-xctr)**2 + (yy-yctr)**2) - - index = r_zone/dx_fine - - vol = (xr**2 - xl**2)*(yr - yl) - - ! weight the zone's data by its size - dens_bin(index) = dens_bin(index) + & - p(i,j,k,dens_comp) * vol - - vel = sqrt(p(i,j,k,xmom_comp)**2 + & - p(i,j,k,ymom_comp)**2) / & - p(i,j,k,dens_comp) - vel_bin(index) = vel_bin(index) + vel * vol - - pres_bin(index) = pres_bin(index) + & - p(i,j,k,pres_comp) * vol - - e_bin(index) = e_bin(index) + & - (p(i,j,k,rhoe_comp) / p(i,j,k,dens_comp)) * vol - - volcount(index) = volcount(index) + vol - - imask(i*r1:(i+1)*r1-1, & - j*r1:(j+1)*r1-1) = 0 - - end if - - enddo - enddo - -end subroutine fextract2d_sph - - -subroutine fextract3d_cyl(lo, hi, p, plo, phi, nc_p, nbins, dens_bin, & - vel_bin, pres_bin, ncount, imask, mask_size, r1, & - dens_comp, xmom_comp, ymom_comp, zmom_comp, pres_comp, dx_fine, dx, & - xctr, yctr) bind(C, name='fextract3d_cyl') - - use amrex_fort_module, only : rt => amrex_real - use amrex_constants_module - - implicit none - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: plo(3), phi(3), nc_p - integer, intent(in), value :: nbins - real(rt), intent(in) :: p(plo(1):phi(1),plo(2):phi(2),plo(3):phi(3),0:nc_p-1) - real(rt), intent(inout) :: dens_bin(0:nbins-1) - real(rt), intent(inout) :: vel_bin(0:nbins-1) - real(rt), intent(inout) :: pres_bin(0:nbins-1) - integer, intent(inout) :: ncount(0:nbins-1) - integer, intent(in), value :: mask_size - integer, intent(inout) :: imask(0:mask_size-1,0:mask_size-1,0:mask_size-1) - integer, intent(in), value :: r1, dens_comp, xmom_comp, ymom_comp, zmom_comp, pres_comp - real(rt), intent(in), value :: dx_fine, xctr, yctr - real(rt), intent(in) :: dx(3) - - integer :: i, j, k, index - real(rt) :: xx, yy, zz, r_zone - - ! loop over all of the zones in the patch. Here, we convert - ! the cell-centered indices at the current level into the - ! corresponding RANGE on the finest level, and test if we've - ! stored data in any of those locations. If we haven't then - ! we store this level's data and mark that range as filled. - do k = lo(3), hi(3) - zz = (k + HALF)*dx(3) - do j = lo(2), hi(2) - yy = (j + HALF)*dx(2) - do i = lo(1), hi(1) - xx = (i + HALF)*dx(1) - - if ( any(imask(i*r1:(i+1)*r1-1, & - j*r1:(j+1)*r1-1, & - k*r1:(k+1)*r1-1) .eq. 1) ) then - - r_zone = sqrt((xx-xctr)**2 + (yy-yctr)**2) - - index = r_zone/dx_fine - - dens_bin(index) = dens_bin(index) + & - p(i,j,k,dens_comp)*r1**3 - - vel_bin(index) = vel_bin(index) + & - (sqrt(p(i,j,k,xmom_comp)**2 + & - p(i,j,k,ymom_comp)**2 + & - p(i,j,k,zmom_comp)**2)/ & - p(i,j,k,dens_comp))*r1**3 - - pres_bin(index) = pres_bin(index) + & - p(i,j,k,pres_comp)*r1**3 - - ncount(index) = ncount(index) + r1**3 - - imask(i*r1:(i+1)*r1-1, & - j*r1:(j+1)*r1-1, & - k*r1:(k+1)*r1-1) = 0 - - end if - - enddo - enddo - enddo - -end subroutine fextract3d_cyl - -subroutine fextract3d_sph(lo, hi, p, plo, phi, nc_p, nbins, dens_bin, & - vel_bin, pres_bin, e_bin, ncount, imask, mask_size, r1,& - dens_comp, xmom_comp, ymom_comp, zmom_comp, pres_comp, rhoe_comp, dx_fine, dx, & - xctr, yctr, zctr) & - bind(C, name='fextract3d_sph') - - use amrex_fort_module, only : rt => amrex_real - use amrex_constants_module - - implicit none - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: plo(3), phi(3), nc_p - integer, intent(in), value :: nbins - real(rt), intent(in) :: p(plo(1):phi(1),plo(2):phi(2),plo(3):phi(3),0:nc_p-1) - real(rt), intent(inout) :: dens_bin(0:nbins-1) - real(rt), intent(inout) :: vel_bin(0:nbins-1) - real(rt), intent(inout) :: pres_bin(0:nbins-1) - real(rt), intent(inout) :: e_bin(0:nbins-1) - integer, intent(inout) :: ncount(0:nbins-1) - integer, intent(in), value :: mask_size - integer, intent(inout) :: imask(0:mask_size-1,0:mask_size-1,0:mask_size-1) - integer, intent(in), value :: r1, dens_comp, xmom_comp, ymom_comp, zmom_comp, pres_comp, rhoe_comp - real(rt), intent(in), value :: dx_fine, xctr, yctr, zctr - real(rt), intent(in) :: dx(3) - - integer :: i, j, k, index - real(rt) :: xx, yy, zz, r_zone, vol, vel - - ! write(*,*) "imask = ", imask - - ! loop over all of the zones in the patch. Here, we convert - ! the cell-centered indices at the current level into the - ! corresponding RANGE on the finest level, and test if we've - ! stored data in any of those locations. If we haven't then - ! we store this level's data and mark that range as filled. - do k = lo(3), hi(3) - zz = (dble(k) + HALF)*dx(3) - do j = lo(2), hi(2) - yy = (dble(j) + HALF)*dx(2) - - do i = lo(1), hi(1) - xx = (dble(i) + HALF)*dx(1) - - if ( any(imask(i*r1:(i+1)*r1-1, & - j*r1:(j+1)*r1-1, & - k*r1:(k+1)*r1-1) .eq. 1) ) then - - r_zone = sqrt((xx-xctr)**2 + (yy-yctr)**2 + (zz-zctr)**2) - - index = r_zone/dx_fine - - ! write(*,*) p(i,1,1,dens_comp) - - ! weight the zone's data by its size - dens_bin(index) = dens_bin(index) + & - p(i,j,k,dens_comp)*r1**3 - - vel_bin(index) = vel_bin(index) + & - (sqrt(p(i,j,k,xmom_comp)**2 + & - p(i,j,k,ymom_comp)**2 + & - p(i,j,k,zmom_comp)**2)/ & - p(i,j,k,dens_comp))*r1**3 - - pres_bin(index) = pres_bin(index) + & - p(i,j,k,pres_comp)*r1**3 - - e_bin(index) = e_bin(index) + & - (p(i,j,k,rhoe_comp)/p(i,j,k,dens_comp))*r1**3 - - ncount(index) = ncount(index) + r1**3 - - imask(i*r1:(i+1)*r1-1, & - j*r1:(j+1)*r1-1, & - k*r1:(k+1)*r1-1) = 0 - - end if - - enddo - enddo - enddo - -end subroutine fextract3d_sph diff --git a/Diagnostics/timestep_limiter/GNUmakefile b/Diagnostics/timestep_limiter/GNUmakefile deleted file mode 100644 index 01103a700f..0000000000 --- a/Diagnostics/timestep_limiter/GNUmakefile +++ /dev/null @@ -1,36 +0,0 @@ -PRECISION = DOUBLE -PROFILE = FALSE -DEBUG = TRUE -DIM = 3 - -COMP = gnu - -USE_MPI = FALSE -USE_OMP = FALSE - -USE_REACT = TRUE - -USE_ACC = FALSE - -# programs to be compiled -ALL: limiter_$(DIM)d.ex - -EOS_DIR := helmholtz - -NETWORK_DIR := aprox21 - -Bpack := ./Make.package -Blocs := . -# EXTERN_SEARCH = . - -CASTRO_HOME := ../.. - -INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Extern/amrdata -include $(AMREX_HOME)/Src/Extern/amrdata/Make.package -vpathdir += $(AMREX_HOME)/Src/Extern/amrdata - -include $(CASTRO_HOME)/Exec/Make.Castro - -limiter_$(DIM)d.ex: $(objForExecs) - @echo Linking $@ ... - $(SILENT) $(PRELINK) $(CXX) $(CPPFLAGS) $(CXXFLAGS) $(LDFLAGS) -o $@ $^ $(libraries) diff --git a/Diagnostics/timestep_limiter/Limiter_F.H b/Diagnostics/timestep_limiter/Limiter_F.H deleted file mode 100644 index 3154c12f22..0000000000 --- a/Diagnostics/timestep_limiter/Limiter_F.H +++ /dev/null @@ -1,41 +0,0 @@ -#ifndef _Limiter_F_H_ -#define _Limiter_F_H_ -#include - - -#ifdef __cplusplus -extern "C" -{ -#endif - -void microphysics_initialize(const int* name, const int* namlen); -void microphysics_finalize(); - -void find_timestep_limiter(const int* lo, const int* hi, - BL_FORT_FAB_ARG_3D(p),const int* nc_p, - const int dens_comp, const int xmom_comp, - const int ymom_comp, const int zmom_comp, - const int pres_comp, const int rhoe_comp, - const int spec_comp, const int time_integration_method, - const amrex::Real* dx, amrex::Real* dt, - amrex::Real* dt_loc); - -void find_timestep_limiter_burning(const int* lo, const int* hi, - BL_FORT_FAB_ARG_3D(p),const int* nc_p, - const int dens_comp, const int temp_comp, - const int rhoe_comp, const int spec_comp, - const amrex::Real* dx, amrex::Real* dt, - amrex::Real* dt_loc); - -void find_timestep_limiter_diffusion(const int* lo, const int* hi, - BL_FORT_FAB_ARG_3D(p),const int* nc_p, - const int dens_comp, const int temp_comp, - const int rhoe_comp, const int spec_comp, - const amrex::Real* dx, amrex::Real* dt, - amrex::Real* dt_loc); - -#ifdef __cplusplus -} -#endif - -#endif diff --git a/Diagnostics/timestep_limiter/Make.package b/Diagnostics/timestep_limiter/Make.package deleted file mode 100644 index 78302a2908..0000000000 --- a/Diagnostics/timestep_limiter/Make.package +++ /dev/null @@ -1,2 +0,0 @@ -ca_F90EXE_sources += limiter_util.F90 -CEXE_headers += Limiter_F.H diff --git a/Diagnostics/timestep_limiter/limiter_util.F90 b/Diagnostics/timestep_limiter/limiter_util.F90 deleted file mode 100644 index 226442167d..0000000000 --- a/Diagnostics/timestep_limiter/limiter_util.F90 +++ /dev/null @@ -1,276 +0,0 @@ -! Process a plotfile to find the location where the timestep is being limited. - -subroutine microphysics_initialize(name, namlen) bind(C, name='microphysics_initialize') - use eos_module, only: eos_init - use network, only: network_init - use actual_rhs_module, only: actual_rhs_init - use amrex_fort_module, only : rt => amrex_real - - integer, intent(in) :: namlen - integer, intent(in) :: name(namlen) - - real(rt) :: small_dens = 1.e-20_rt, small_temp = 1.e-20_rt - - call runtime_init(name, namlen) - - call eos_init(small_dens=small_dens, small_temp=small_temp) - - call network_init() - call actual_rhs_init() - -end subroutine microphysics_initialize - -subroutine microphysics_finalize() bind(C, name='microphysics_finalize') - - use eos_module, only: eos_finalize - use network, only: network_finalize - - call eos_finalize() - call network_finalize() - -end subroutine microphysics_finalize - -subroutine find_timestep_limiter(lo, hi, state, slo, shi, nc_s, & - dens_comp, xmom_comp, ymom_comp, zmom_comp, pres_comp, rhoe_comp, spec_comp, time_integration_method, dx, dt, dt_loc) & - bind(C, name='find_timestep_limiter') - - use amrex_fort_module, only : rt => amrex_real, amrex_min - use amrex_constants_module - use eos_module, only: eos - use eos_type_module, only: eos_t, eos_input_re - use network, only: nspec - - implicit none - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: slo(3), shi(3), nc_s - real(rt), intent(in) :: state(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3),0:nc_s-1) - integer, intent(in), value :: dens_comp, xmom_comp, ymom_comp, zmom_comp, pres_comp, rhoe_comp, spec_comp, time_integration_method - real(rt), intent(in) :: dx(3) - real(rt), intent(inout) :: dt - real(rt), intent(inout) :: dt_loc(3) - - integer :: i, j, k - real(rt) :: x, y, z, cs, dt1, dt2, dt3, dt_tmp, ux, uy, uz - - type (eos_t) :: eos_state - - do k = lo(3), hi(3) - z = (dble(k) + HALF)*dx(3) - do j = lo(2), hi(2) - y = (dble(j) + HALF)*dx(2) - - do i = lo(1), hi(1) - x = (dble(i) + HALF)*dx(1) - - eos_state % rho = state(i,j,k,dens_comp) - eos_state % e = state(i,j,k,rhoe_comp) / eos_state % rho - eos_state % xn = state(i,j,k,spec_comp:spec_comp+nspec-1) - - call eos(eos_input_re, eos_state) - - cs = eos_state % cs - - ! calculate CFL timestep using the velocity - ux = state(i,j,k,xmom_comp) / state(i,j,k,dens_comp) - - dt1 = dx(1) / (cs + abs(ux)) - dt2 = dt1 - dt3 = dt1 -#if AMREX_SPACEDIM >= 2 - uy = state(i,j,k,ymom_comp) / state(i,j,k,dens_comp) - dt2 = dx(2) / (cs + abs(uy)) -#endif -#if AMREX_SPACEDIM == 3 - uz = state(i,j,k,zmom_comp) / state(i,j,k,dens_comp) - dt3 = dx(3) / (cs + abs(uz)) -#endif - - if (time_integration_method == 0) then - if (min(dt1, dt2, dt3) < dt) then - dt = min(dt1, dt2, dt3) - dt_loc(1:3) = (/ x, y, z /) - endif - else - dt_tmp = ONE / dt1 -#if AMREX_SPACEDIM >= 2 - dt_tmp = dt_tmp + ONE / dt2 -#endif -#if AMREX_SPACEDIM == 3 - dt_tmp = dt_tmp + ONE / dt3 -#endif - if (ONE/dt_tmp < dt) then - dt = ONE/dt_tmp - dt_loc(1:3) = (/ x, y, z /) - endif - - endif - enddo - enddo - enddo - -end subroutine find_timestep_limiter - - -subroutine find_timestep_limiter_burning(lo, hi, state, slo, shi, nc_s, & - dens_comp, temp_comp, rhoe_comp, spec_comp, dx, dt, dt_loc) & - bind(C, name='find_timestep_limiter_burning') - - use amrex_fort_module, only : rt => amrex_real, amrex_min - use amrex_constants_module - use actual_rhs_module, only: actual_rhs - use eos_module, only: eos - use eos_type_module, only: eos_t, eos_input_rt - use burner_module, only: ok_to_burn - use burn_type_module, only : burn_t, net_ienuc, burn_to_eos, eos_to_burn - use network, only: nspec, aion - use meth_params_module, only : dtnuc_e, dtnuc_X, dtnuc_X_threshold - use extern_probin_module, only: small_x - - implicit none - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: slo(3), shi(3), nc_s - real(rt), intent(in) :: state(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3),0:nc_s-1) - integer, intent(in), value :: dens_comp, temp_comp, rhoe_comp, spec_comp - real(rt), intent(in) :: dx(3) - real(rt), intent(inout) :: dt - real(rt), intent(inout) :: dt_loc(3) - - integer :: i, j, k, n - real(rt) :: x, y, z, e, Xn(nspec), dedt, dXdt(nspec) - - type (eos_t) :: eos_state - - type (burn_t) :: state_old, state_new - - ! Set a floor on the minimum size of a derivative. This floor - ! is small enough such that it will result in no timestep limiting. - - real(rt), parameter :: derivative_floor = 1.e-50_rt - - do k = lo(3), hi(3) - z = (dble(k) + HALF)*dx(3) - do j = lo(2), hi(2) - y = (dble(j) + HALF)*dx(2) - do i = lo(1), hi(1) - x = (dble(i) + HALF)*dx(1) - - state_old % rho = state(i,j,k,dens_comp) - state_old % T = state(i,j,k,temp_comp) - state_old % e = state(i,j,k,rhoe_comp) / state(i,j,k,dens_comp) - state_old % xn = state(i,j,k,spec_comp:spec_comp+nspec-1) - - state_new % rho = state(i,j,k,dens_comp) - state_new % T = state(i,j,k,temp_comp) - state_new % e = state(i,j,k,rhoe_comp) / state(i,j,k,dens_comp) - state_new % xn = state(i,j,k,spec_comp:spec_comp+nspec-1) - - ! if (.not. ok_to_burn(state_new)) cycle - - e = state_new % e - Xn = max(state_new % xn, small_x) - - call burn_to_eos(state_new, eos_state) - call eos(eos_input_rt, eos_state) - call eos_to_burn(eos_state, state_new) - - state_new % dx = minval(dx(1:AMREX_SPACEDIM)) - - call actual_rhs(state_new) - - dedt = state_new % ydot(net_ienuc) - dXdt = state_new % ydot(1:nspec) * aion - - dedt = max(abs(dedt), derivative_floor) - - do n = 1, nspec - if (Xn(n) .ge. dtnuc_X_threshold) then - dXdt(n) = max(abs(dXdt(n)), derivative_floor) - else - dXdt(n) = derivative_floor - endif - enddo - - if (dtnuc_e * e / dedt < dt) then - dt = dtnuc_e * e / dedt - dt_loc = (/ x, y, z /) - else if (dtnuc_X * minval(X / dXdt) < dt) then - dt = dtnuc_X * minval(X / dXdt) - dt_loc = (/ x, y, z /) - endif - - enddo - enddo - enddo - -end subroutine find_timestep_limiter_burning - -#ifdef DIFFUSION -subroutine find_timestep_limiter_diffusion(lo, hi, state, slo, shi, nc_s, & - dens_comp, temp_comp, rhoe_comp, spec_comp, dx, dt, dt_loc) & - bind(C, name='find_timestep_limiter_diffusion') - - use amrex_fort_module, only : rt => amrex_real, amrex_min - use amrex_constants_module - use eos_module, only: eos - use eos_type_module, only: eos_t, eos_input_re - use conductivity_module, only: conductivity - use meth_params_module, only : diffuse_cutoff_density - - implicit none - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: slo(3), shi(3), nc_s - real(rt), intent(in) :: state(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3),0:nc_s-1) - integer, intent(in), value :: dens_comp, temp_comp, rhoe_comp, spec_comp - real(rt), intent(in) :: dx(3) - real(rt), intent(inout) :: dt - real(rt), intent(inout) :: dt_loc(3) - - integer :: i, j, k - real(rt) :: x, y, z, dt1, dt2, dt3, D - - type (eos_t) :: eos_state - - do k = lo(3), hi(3) - z = (dble(k) + HALF)*dx(3) - do j = lo(2), hi(2) - y = (dble(j) + HALF)*dx(2) - do i = lo(1), hi(1) - x = (dble(i) + HALF)*dx(1) - - if (state(i,j,k,dens_comp) > diffuse_cutoff_density) then - - eos_state % rho = state(i,j,k,dens_comp) - eos_state % T = state(i,j,k,temp_comp) - eos_state % e = state(i,j,k,rhoe_comp) / state(i,j,k,dens_comp) - state_old % xn = state(i,j,k,spec_comp:spec_comp+nspec-1) - - call eos(eos_input_re, eos_state) - - call conductivity(eos_state) - - D = eos_state % conductivity / (eos_state % rho * eos_state % cv) - - dt1 = HALF * dx(1)**2 / D - dt2 = dt1 - dt3 = dt1 - -#if AMREX_SPACEDIM >= 2 - dt2 = HALF * dx(2)**2 / D -#endif -#if AMREX_SPACEDIM == 3 - dt3 = HALF * dx(3)**2 / D -#endif - if (min(dt1, dt2, dt3) < dt) then - dt = min(dt1, dt2, dt3) - dt_loc = (/ x, y, z /) - endif - endif - enddo - enddo - enddo - -end subroutine find_timestep_limiter_diffusion -#endif diff --git a/Diagnostics/timestep_limiter/main.cpp b/Diagnostics/timestep_limiter/main.cpp deleted file mode 100644 index a10642975d..0000000000 --- a/Diagnostics/timestep_limiter/main.cpp +++ /dev/null @@ -1,271 +0,0 @@ -// -// Process a plotfile to find the location where the timestep is being limited -// -#include -#include -// #include -#include -#include -#include -#include -#include -#include - -using namespace amrex; - -std::string inputs_name = ""; - -// -// Prototypes -// -void GetInputArgs (const int argc, char** argv, - string& pltfile); - -void ProcessJobInfo(string job_info_file, string inputs_file_name); - - -int main(int argc, char* argv[]) -{ - - int dummy = 0; - - amrex::Initialize(dummy, argv); - - // timer for profiling - BL_PROFILE_VAR("main()", pmain); - - // Tell AMReX to be quiet - amrex::system::verbose = 0; - - // Input arguments - string pltfile; - - GetInputArgs (argc, argv, pltfile); - - string job_info = pltfile + "/job_info"; - string inputs_file = "inputs.txt"; - - ProcessJobInfo(job_info, inputs_file); - - // this will process the input parameters from the job_info file - amrex::ParmParse::Initialize(0,0,inputs_file.c_str()); - - // Read in the input values to Fortran. - int NUM_GROW; - - ca_get_method_params(&NUM_GROW); - ca_set_castro_method_params(); - - // Start dataservices - DataServices::SetBatchMode(); - - // Define the type of file - Amrvis::FileType fileType(Amrvis::NEWPLT); - DataServices dataServices (pltfile, fileType); - - if (!dataServices.AmrDataOk()) - DataServices::Dispatch(DataServices::ExitRequest, NULL); - - // initialize microphysics stuff - auto probin_name = "probin"; - auto probin_file_length = 6; - Vector probin_file(probin_file_length); - for (auto i = 0; i < probin_file_length; i++) { - probin_file[i] = probin_name[i]; - } - - microphysics_initialize(probin_file.dataPtr(), &probin_file_length); - - // get data from plot file - AmrData& data = dataServices.AmrDataRef(); - - // get variable names - const Vector& varNames = data.PlotVarNames(); - - int dens_comp, xmom_comp, ymom_comp, zmom_comp, pres_comp, rhoe_comp, spec_comp, temp_comp; - int time_integration_method = 3; - - // find variable indices - dens_comp = data.StateNumber("density"); - xmom_comp = data.StateNumber("xmom"); -#if (AMREX_SPACEDIM >= 2) - ymom_comp = data.StateNumber("ymom"); -#endif -#if (AMREX_SPACEDIM == 3) - zmom_comp = data.StateNumber("zmom"); -#endif - pres_comp = data.StateNumber("pressure"); - rhoe_comp = data.StateNumber("rho_e"); - temp_comp = data.StateNumber("Temp"); - - if (dens_comp < 0 || xmom_comp < 0 || pres_comp < 0 || rhoe_comp < 0) - Abort("ERROR: variable(s) not found"); - -#if (AMREX_SPACEDIM == 3) - if (ymom_comp < 0 || zmom_comp < 0) - Abort("ERROR: variable(s) not found"); -#endif - - // we're going to find the spec comp by looking for the first variable name - // that begins with 'X' - std::string first_spec_name = ""; - for (auto &it : varNames) { - if (it[0] == 'X') { - first_spec_name = it; - break; - } - } - - if (first_spec_name == "") - Abort("ERROR: no species were found"); - - spec_comp = data.StateNumber(first_spec_name); - - // fill a multifab with the data - Vector fill_comps(data.NComp()); - for (auto i = 0; i < data.NComp(); i++) - fill_comps[i] = i; - - Vector dt_loc = {0.,0.,0.}; - Vector burning_dt_loc = {0.,0.,0.}; - Vector diffusion_dt_loc = {0.,0.,0.}; - Real dt = 1.e99; - Real burning_dt = 1.e99; - Real diffusion_dt = 1.e99; - - // loop over the data - for (int lev=0; lev <= data.FinestLevel(); lev++) { - - Vector level_dx = data.DxLevel()[lev]; - - const BoxArray& ba = data.boxArray(lev); - const DistributionMapping& dm = data.DistributionMap(lev); - - MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); - - for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); - - // we only want to load the data a Fab at the time as otherwise - // we can experience memory issues. - // Therefore, we'll iterate over the variables using FillVar - // to fill a temporary fab, then copy this over the to MultiFab - Array4 lev_data_arr = lev_data_mf.array(mfi); - FArrayBox tmp_fab(bx); - - for (auto n = 0; n < data.NComp(); n++) { - data.FillVar(&tmp_fab, bx, lev, varNames[n], ParallelDescriptor::IOProcessorNumber()); - - Array4 tmp_arr = tmp_fab.array(); - - AMREX_PARALLEL_FOR_3D(bx, i, j, k, { - lev_data_arr(i,j,k,n) = tmp_arr(i,j,k); - }); - } - - find_timestep_limiter(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - dens_comp, xmom_comp, ymom_comp, zmom_comp, pres_comp, rhoe_comp, spec_comp, - time_integration_method, - level_dx.dataPtr(), &dt, dt_loc.dataPtr()); - - find_timestep_limiter_burning(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - dens_comp, temp_comp, rhoe_comp, spec_comp, - level_dx.dataPtr(), &burning_dt, burning_dt_loc.dataPtr()); -#ifdef DIFFUSION - find_timestep_limiter_diffusion(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - dens_comp, temp_comp, rhoe_comp, spec_comp, - level_dx.dataPtr(), &diffusion_dt, diffusion_dt_loc.dataPtr()); -#endif - } - } - Print() << std::endl; - - Print() << "dt = " << dt << " at location"; - for (auto i = 0; i < AMREX_SPACEDIM; i++) { - Print() << ' ' << dt_loc[i]; - } - Print() << std::endl; - - Print() << "burning_dt = " << burning_dt << " at location"; - for (auto i = 0; i < AMREX_SPACEDIM; i++) { - Print() << ' ' << burning_dt_loc[i]; - } - Print() << std::endl; - -#ifdef DIFFUSION - Print() << "diffusion_dt = " << diffusion_dt << " at location"; - for (auto i = 0; i < AMREX_SPACEDIM; i++) { - Print() << ' ' << diffusion_dt_loc[i]; - } - Print() << std::endl; -#endif - - Print() << std::endl; - - // finalize microphysics stuff - microphysics_finalize(); - - // destroy timer for profiling - BL_PROFILE_VAR_STOP(pmain); - - amrex::Finalize(); -} - -// -// Parse command line arguments -// -void GetInputArgs ( const int argc, char** argv, - string& pltfile) -{ - pltfile = argv[1]; - - if (pltfile.empty()) - { - Abort("Missing input file"); - } - - Print() << "Finding limiting timestep in plotfile = \"" << pltfile << "\"" << std::endl; -} - -// -// Reads in a job_info file, extracts the inputs parameters and saves -// them to a new file -// -void ProcessJobInfo(string job_info_file, string inputs_file_name) -{ - std::ifstream job_info (job_info_file); - std::ofstream inputs_file (inputs_file_name); - - std::regex inputs_rgx("Inputs File Parameters"); - - if (job_info.is_open() && inputs_file.is_open()) { - string line; - bool found_inputs = false; - while (getline(job_info, line)) { - if (found_inputs) { - // if the last character is a space, then there is no variable - if (line.back() != ' ') { - if (line.front() == '[') { - inputs_file << line.substr(3) << std::endl; - } else { - inputs_file << line << std::endl; - } - } - } else { - if (std::regex_search(line, inputs_rgx)) { - // found a match! - found_inputs = true; - - // skip a line - getline(job_info, line); - } - } - } - - job_info.close(); - inputs_file.close(); - } -} \ No newline at end of file diff --git a/Diagnostics/timestep_limiter/probin b/Diagnostics/timestep_limiter/probin deleted file mode 100644 index c3541df6af..0000000000 --- a/Diagnostics/timestep_limiter/probin +++ /dev/null @@ -1,3 +0,0 @@ -&extern - -/ diff --git a/Docs/source/FlowChart.rst b/Docs/source/FlowChart.rst index 90c65f3969..d91c3718c9 100644 --- a/Docs/source/FlowChart.rst +++ b/Docs/source/FlowChart.rst @@ -322,7 +322,8 @@ In the code, the objective is to evolve the state from the old time, A. In ``initialize_do_advance()``, create ``Sborder``, initialized from ``S_old`` - B. Check for NaNs in the initial state, ``S_old``. + B. Call ``clean_state()`` to make sure the thermodynamics are in sync, in particular, + compute the temperature. #. *React* :math:`\Delta t/2` [``strang_react_first_half()`` ] @@ -391,7 +392,12 @@ In the code, the objective is to evolve the state from the old time, rather than computing it from the Riemann problem. This source is computed here for the internal energy equation. - D. [``DIFFUSION``] diffusion : thermal diffusion can be + D. geometry source: this is applied only for 2-d axisymmetric data + and captures the geometric term arising from applying the + cylindrical divergence in :math:`\nabla \cdot (\rho \Ub \Ub)` in + the momentum equation. See :cite:`bernard-champmartin_eulerian_2012`. + + E. [``DIFFUSION``] diffusion : thermal diffusion can be added in an explicit formulation. Second-order accuracy is achieved by averaging the time-level :math:`n` and :math:`n+1` terms, using the same predictor-corrector strategy described here. @@ -402,10 +408,10 @@ In the code, the objective is to evolve the state from the old time, timestep constraint, since the treatment is explicit. See Chapter :ref:`ch:diffusion` for more details. - E. [``HYBRID_MOMENTUM``] angular momentum + F. [``HYBRID_MOMENTUM``] angular momentum - F. [``GRAVITY``] gravity: + G. [``GRAVITY``] gravity: For full Poisson gravity, we solve for for gravity using: @@ -420,7 +426,7 @@ In the code, the objective is to evolve the state from the old time, solver are given in Chapter :ref:`ch:gravity`. - G. [``ROTATION``] rotation + H. [``ROTATION``] rotation We compute the rotational potential (for use in the energy update) and the rotational acceleration (for use in the momentum diff --git a/Docs/source/refs.bib b/Docs/source/refs.bib index 45b1eae848..6d2f07cd6e 100644 --- a/Docs/source/refs.bib +++ b/Docs/source/refs.bib @@ -1063,3 +1063,35 @@ @article{harpole2021dynamics subject = {X-ray bursts}, url = {https://ui.adsabs.harvard.edu/abs/2021arXiv210200051H/abstract} } + +@article{bernard-champmartin_eulerian_2012, + title = {An {Eulerian} finite volume solver for multi-material fluid flows with cylindrical symmetry.}, + url = {https://hal.archives-ouvertes.fr/hal-00797200}, + doi = {10.1016/j.compfluid.2012.09.014}, + abstract = {In this paper, we adapt a pre-existing 2D + cartesian cell centered finite volume solver to + treat the compressible 3D Euler equations with + cylindrical symmetry. We the n extend it to + multi-material flows. Assuming cylindrical symmetry + with respect to the z axis (i.e. all the functions + do not depend explicitly on the angular variable + \${\textbackslash} theta\$), we obtain a set of five + conservation laws with source terms that can be + decoupled in two systems solved on a 2D orthogonal + mesh in which a cell as a torus geometry. A specific + upwinding treatment of the source term is required + and implemented for the stationary case. Test cases + will be presented for vanishing and non-vanishing + azimuthal velocity \$u\_\{{\tex + tbackslash}theta\}\$.}, + urldate = {2020-03-11}, + journal = {Computers and Fluids}, + author = {Bernard-Champmartin, Aude and Braeunig, Jean-Philippe and Ghidaglia, Jean-Michel}, + month = sep, + year = {2012}, + note = {Publisher: Elsevier}, + keywords = {Axisymmetric compressible Euler equations, Finite volumes, Multi-material flows}, + pages = {7}, + file = {HAL PDF Full Text:/home/zingale/Zotero/storage/FZANNPPD/Bernard-Champmartin et al. - 2012 - An Eulerian finite volume solver for multi-materia.pdf:application/pdf}, +} + diff --git a/Exec/Make.auto_source b/Exec/Make.auto_source index cd3bb87eef..88248e98dd 100644 --- a/Exec/Make.auto_source +++ b/Exec/Make.auto_source @@ -16,11 +16,8 @@ CEXE_headers += extern_parameters.H # for dependency resolution AUTO_BUILD_SOURCES += $(CASTRO_AUTO_SOURCE_DIR)/extern_parameters.H -# these are for finding runtime parameters -EXTERN_SEARCH += $(EXTERN_CORE) -EXTERN_SEARCH += $(MICROPHYSICS_HOME)/networks/ - -EXTERN_PARAMETERS := $(foreach dir, $(EXTERN_SEARCH),$(realpath $(wildcard $(dir)/_parameters))) +# this is for finding runtime parameters +EXTERN_PARAMETERS := $(foreach dir, $(INCLUDE_LOCATIONS),$(realpath $(wildcard $(dir)/_parameters))) $(CASTRO_AUTO_SOURCE_DIR)/extern_parameters.cpp: $(CASTRO_AUTO_SOURCE_DIR)/extern_parameters.H diff --git a/Exec/gravity_tests/hydrostatic_adjust/Problem_Derive.cpp b/Exec/gravity_tests/hydrostatic_adjust/Problem_Derive.cpp index bbbd46b0e4..0c9058f24d 100644 --- a/Exec/gravity_tests/hydrostatic_adjust/Problem_Derive.cpp +++ b/Exec/gravity_tests/hydrostatic_adjust/Problem_Derive.cpp @@ -55,7 +55,7 @@ void ca_derpi(const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, for (int n = 0; n < NumSpec; n++) { eos_state.xn[n] = dat(i,j,k,UFS+n) / dat(i,j,k,URHO); } -#if NAUX > 0 +#if NAUX_NET > 0 for (int n = 0; n < NumAux; n++) { eos_state.aux[n] = dat(i,j,k,UFX+n) / dat(i,j,k,URHO); } @@ -132,7 +132,7 @@ void ca_derpioverp0(const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, for (int n = 0; n < NumSpec; n++) { eos_state.xn[n] = dat(i,j,k,UFS+n) / dat(i,j,k,URHO); } -#if NAUX > 0 +#if NAUX_NET > 0 for (int n = 0; n < NumAux; n++) { eos_state.aux[n] = dat(i,j,k,UFX+n) / dat(i,j,k,URHO); } diff --git a/Exec/hydro_tests/gamma_law_bubble/Problem_Derive.cpp b/Exec/hydro_tests/gamma_law_bubble/Problem_Derive.cpp index 5c9e44e3a7..7f88ab464e 100644 --- a/Exec/hydro_tests/gamma_law_bubble/Problem_Derive.cpp +++ b/Exec/hydro_tests/gamma_law_bubble/Problem_Derive.cpp @@ -55,7 +55,7 @@ void ca_derpi(const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, for (int n = 0; n < NumSpec; n++) { eos_state.xn[n] = dat(i,j,k,UFS+n) / dat(i,j,k,URHO); } -#if NAUX > 0 +#if NAUX_NET > 0 for (int n = 0; n < NumAux; n++) { eos_state.aux[n] = dat(i,j,k,UFX+n) / dat(i,j,k,URHO); } @@ -120,7 +120,7 @@ void ca_derpioverp0(const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, for (int n = 0; n < NumSpec; n++) { eos_state.xn[n] = dat(i,j,k,UFS+n) * rhoInv; } -#if NAUX > 0 +#if NAUX_NET > 0 for (int n = 0; n < NumAux; n++) { eos_state.aux[n] = dat(i,j,k,UFX+n) * rhoInv; } diff --git a/Exec/reacting_tests/nse_test/GNUmakefile b/Exec/reacting_tests/nse_test/GNUmakefile index 36c1de76a2..0c9593102e 100644 --- a/Exec/reacting_tests/nse_test/GNUmakefile +++ b/Exec/reacting_tests/nse_test/GNUmakefile @@ -19,7 +19,7 @@ EOS_DIR := helmholtz # This sets the EOS directory in $(MICROPHYSICS_HOME)/networks NETWORK_DIR := aprox19 -USE_NSE := TRUE +USE_NSE_TABLE := TRUE INTEGRATOR_DIR := VODE diff --git a/Exec/reacting_tests/nse_test/_prob_params b/Exec/reacting_tests/nse_test/_prob_params index 1735b7830c..4ca9c39649 100644 --- a/Exec/reacting_tests/nse_test/_prob_params +++ b/Exec/reacting_tests/nse_test/_prob_params @@ -1,9 +1,20 @@ - +# ambient density rho0 real 1.4_rt y +# ambient temperature T0 real 1.0_rt y +# amplitude of temperature perturbation dT_fact real 1.5_rt y +# scale of the temperature perturbation L_pert real 0.5_rt y +# ambient x-velocity +u0 real 0.0_rt y + +# ambient y-velocity field +v0 real 0.0_rt y + +# ambient z-velocity field +w0 real 0.0_rt y \ No newline at end of file diff --git a/Exec/reacting_tests/nse_test/analysis/slice_multi.py b/Exec/reacting_tests/nse_test/analysis/slice_multi.py new file mode 100755 index 0000000000..e567c8c4e9 --- /dev/null +++ b/Exec/reacting_tests/nse_test/analysis/slice_multi.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python3 + +import os +import sys +import yt +import matplotlib.pyplot as plt + +from mpl_toolkits.axes_grid1 import ImageGrid +from matplotlib.ticker import ScalarFormatter + +# assume that our data is in CGS +from yt.units import cm + +plotfile = sys.argv[1] + +# slice plot of temperature +ds = yt.load(plotfile, hint="castro") + +xctr = 0.5*(ds.domain_left_edge[0] + ds.domain_right_edge[0]) +L_x = ds.domain_right_edge[0] - ds.domain_left_edge[0] + +yctr = 0.5*(ds.domain_left_edge[1] + ds.domain_right_edge[1]) +L_y = ds.domain_right_edge[1] - ds.domain_left_edge[1] + + +fig = plt.figure() + +nrows = 2 +ncols = 3 + +grid = ImageGrid(fig, 111, nrows_ncols=(nrows, ncols), axes_pad=0.75, cbar_pad="2%", + label_mode="L", cbar_mode="each") + + +fields = ["Ye", "abar", "enuc", "Temp", "MachNumber"] + +for i in range(nrows * ncols): + + if i >= len(fields): + grid[i].remove() + grid.cbar_axes[i].remove() + continue + else: + f = fields[i] + + sp = yt.SlicePlot(ds, "z", f, center=[xctr, yctr, 0.0*cm], + width=[L_x, L_y, 0.0*cm], fontsize="12") + + sp.set_buff_size((2000,2000)) + + if f == "Ye": + sp.set_zlim(f, 0.48, 0.5) + sp.set_log(f, False) + sp.set_cmap(f, "magma_r") + elif f == "abar": + sp.set_log(f, False) + sp.set_cmap(f, "viridis") + elif f == "enuc": + sp.set_log(f, True, linthresh=1.e10) + sp.set_zlim(f, -1.e18, 1.e18) + sp.set_cmap(f, "bwr") + elif f == "MachNumber": + sp.set_zlim(f, 1.e-5, 0.1) + sp.set_cmap(f, "plasma") + elif f == "magvel": + sp.set_zlim(f, 100.0, 2.e7) + sp.set_cmap(f, "viridis") + elif f == "Temp": + sp.set_cmap(f, "magma") + sp.set_zlim(f, 1.e8, 6.e9) + + sp.set_axes_unit("cm") + + plot = sp.plots[f] + plot.figure = fig + plot.axes = grid[i].axes + plot.cax = grid.cbar_axes[i] + if i < len(fields)-1: + grid[i].axes.xaxis.offsetText.set_visible(False) + #fmt = grid.cbar_axes[i].yaxis.set_major_formatter(ScalarFormatter(useMathText=True)) + sp._setup_plots() + + + +fig.set_size_inches(12.0, 9.0) +plt.tight_layout() +plt.savefig("{}_slice.png".format(os.path.basename(plotfile))) + diff --git a/Exec/reacting_tests/nse_test/convergence_simplified_sdc.sh b/Exec/reacting_tests/nse_test/convergence_simplified_sdc.sh index 1110e21f6b..87067a7896 100755 --- a/Exec/reacting_tests/nse_test/convergence_simplified_sdc.sh +++ b/Exec/reacting_tests/nse_test/convergence_simplified_sdc.sh @@ -10,12 +10,12 @@ castro.time_integration_method=3 mpiexec -n 8 ${EXEC} inputs.64 ${RUNPARAMS} >& /dev/null mpiexec -n 8 ${EXEC} inputs.128 ${RUNPARAMS} >& /dev/null -mpiexec -n 8 ${EXEC} inputs.256 ${RUNPARAMS} >& /dev/null +mpiexec -n 16 ${EXEC} inputs.256 ${RUNPARAMS} >& /dev/null -RichardsonConvergenceTest2d.gnu.ex coarFile=nse_test_64_plt00080 mediFile=nse_test_128_plt00160 fineFile=nse_test_256_plt00320 >& nse_convergence_simple_sdc_lo.out +RichardsonConvergenceTest2d.gnu.ex coarFile=nse_test_64_plt00125 mediFile=nse_test_128_plt00250 fineFile=nse_test_256_plt00500 >& nse_convergence_simple_sdc_lo.out -mpiexec -n 8 ${EXEC} inputs.512 ${RUNPARAMS} >& /dev/null +mpiexec -n 16 ${EXEC} inputs.512 ${RUNPARAMS} >& /dev/null -RichardsonConvergenceTest2d.gnu.ex coarFile=nse_test_128_plt00160 mediFile=nse_test_256_plt00320 fineFile=nse_test_512_plt00640 >& nse_convergence_simple_sdc_hi.out +RichardsonConvergenceTest2d.gnu.ex coarFile=nse_test_128_plt00250 mediFile=nse_test_256_plt00500 fineFile=nse_test_512_plt01000 >& nse_convergence_simple_sdc_hi.out diff --git a/Exec/reacting_tests/nse_test/convergence_simplified_sdc_w_vel.sh b/Exec/reacting_tests/nse_test/convergence_simplified_sdc_w_vel.sh new file mode 100755 index 0000000000..8646409ce5 --- /dev/null +++ b/Exec/reacting_tests/nse_test/convergence_simplified_sdc_w_vel.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +# echo the commands +set -x + +EXEC=./Castro2d.gnu.MPI.SMPLSDC.ex +RUNPARAMS=" +castro.time_integration_method=3 +problem.u0=1.e8 +problem.v0=1.e8 +" + +mpiexec -n 8 ${EXEC} inputs.64 ${RUNPARAMS} >& /dev/null +mpiexec -n 8 ${EXEC} inputs.128 ${RUNPARAMS} >& /dev/null +mpiexec -n 8 ${EXEC} inputs.256 ${RUNPARAMS} >& /dev/null + +RichardsonConvergenceTest2d.gnu.ex coarFile=nse_test_64_plt00080 mediFile=nse_test_128_plt00160 fineFile=nse_test_256_plt00320 >& nse_convergence_simple_sdc_lo.out + +mpiexec -n 8 ${EXEC} inputs.512 ${RUNPARAMS} >& /dev/null + +RichardsonConvergenceTest2d.gnu.ex coarFile=nse_test_128_plt00160 mediFile=nse_test_256_plt00320 fineFile=nse_test_512_plt00640 >& nse_convergence_simple_sdc_hi.out + + diff --git a/Exec/reacting_tests/nse_test/convergence_strang.sh b/Exec/reacting_tests/nse_test/convergence_strang.sh index 8829713661..2dae53146f 100755 --- a/Exec/reacting_tests/nse_test/convergence_strang.sh +++ b/Exec/reacting_tests/nse_test/convergence_strang.sh @@ -3,16 +3,17 @@ # echo the commands set -x -DIM=1 +DIM=2 EXEC=./Castro${DIM}d.gnu.MPI.ex -mpiexec -n 4 ${EXEC} inputs.64 >& /dev/null -mpiexec -n 4 ${EXEC} inputs.128 >& /dev/null -mpiexec -n 4 ${EXEC} inputs.256 >& /dev/null -mpiexec -n 4 ${EXEC} inputs.512 >& /dev/null +mpiexec -n 8 ${EXEC} inputs.64 >& /dev/null +mpiexec -n 8 ${EXEC} inputs.128 >& /dev/null +mpiexec -n 8 ${EXEC} inputs.256 >& /dev/null -RichardsonConvergenceTest${DIM}d.gnu.ex coarFile=nse_test_64_plt00080 mediFile=nse_test_128_plt00160 fineFile=nse_test_256_plt00320 >& nse_convergence_strang_lo.out +RichardsonConvergenceTest${DIM}d.gnu.ex coarFile=nse_test_64_plt00125 mediFile=nse_test_128_plt00250 fineFile=nse_test_256_plt00500 >& nse_convergence_strang_lo.out -RichardsonConvergenceTest${DIM}d.gnu.ex coarFile=nse_test_128_plt00160 mediFile=nse_test_256_plt00320 fineFile=nse_test_512_plt00640 >& nse_convergence_strang_hi.out +mpiexec -n 8 ${EXEC} inputs.512 >& /dev/null + +RichardsonConvergenceTest${DIM}d.gnu.ex coarFile=nse_test_128_plt00250 mediFile=nse_test_256_plt00500 fineFile=nse_test_512_plt01000 >& nse_convergence_strang_hi.out diff --git a/Exec/reacting_tests/nse_test/convergence_strang_w_vel.sh b/Exec/reacting_tests/nse_test/convergence_strang_w_vel.sh new file mode 100755 index 0000000000..7ac47f638f --- /dev/null +++ b/Exec/reacting_tests/nse_test/convergence_strang_w_vel.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +# echo the commands +set -x + +DIM=2 +EXEC=./Castro${DIM}d.gnu.MPI.ex + +RUNPARAMS=" +problem.u0=1.e8 +problem.v0=1.e8 +" +mpiexec -n 4 ${EXEC} inputs.64 ${RUNPARAMS} >& /dev/null +mpiexec -n 4 ${EXEC} inputs.128 ${RUNPARAMS} >& /dev/null +mpiexec -n 4 ${EXEC} inputs.256 ${RUNPARAMS} >& /dev/null + +RichardsonConvergenceTest${DIM}d.gnu.ex coarFile=nse_test_64_plt00080 mediFile=nse_test_128_plt00160 fineFile=nse_test_256_plt00320 >& nse_convergence_strang_lo.out + +mpiexec -n 4 ${EXEC} inputs.512 ${RUNPARAMS} >& /dev/null + +RichardsonConvergenceTest${DIM}d.gnu.ex coarFile=nse_test_128_plt00160 mediFile=nse_test_256_plt00320 fineFile=nse_test_512_plt00640 >& nse_convergence_strang_hi.out + + diff --git a/Exec/reacting_tests/nse_test/create_pretty_tables.py b/Exec/reacting_tests/nse_test/create_pretty_tables.py index b94e5c27a7..79e39896b9 100644 --- a/Exec/reacting_tests/nse_test/create_pretty_tables.py +++ b/Exec/reacting_tests/nse_test/create_pretty_tables.py @@ -123,7 +123,8 @@ def read_convergence(file_lo, file_hi): "rho_Fe56": r"$\rho X(\isotm{Fe}{56})$", "rho_Ye": r"$\rho Y_e$", "rho_abar": r"$\rho \bar{A}$", - "rho_bea": r"$\rho (B/A)$" + "rho_bea": r"$\rho (B/A)$", + "rho_enuc": r"$\rho \dot{S}$" } # sdc4 @@ -134,5 +135,9 @@ def read_convergence(file_lo, file_hi): for v in sdc4.data: if v.name in good_vars.keys(): - print(v.get_table_line(pretty_name=good_vars[v.name], simple=args.simple)) + if args.simple: + name = v.name + else: + name = good_vars[v.name] + print(v.get_table_line(pretty_name=name, simple=args.simple)) diff --git a/Exec/reacting_tests/nse_test/inputs.128 b/Exec/reacting_tests/nse_test/inputs.128 index e9ef58a589..eec07b0c24 100644 --- a/Exec/reacting_tests/nse_test/inputs.128 +++ b/Exec/reacting_tests/nse_test/inputs.128 @@ -1,12 +1,12 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- max_step = 15000 -stop_time = 0.08 +stop_time = 0.025 # PROBLEM SIZE & GEOMETRY geometry.is_periodic = 1 1 1 geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical geometry.prob_lo = 0 0 0 -geometry.prob_hi = 1.e8 1.e8 1.e8 +geometry.prob_hi = 2.e7 2.e7 2.e7 amr.n_cell = 128 128 128 @@ -35,7 +35,7 @@ castro.small_temp = 1.e7 castro.cfl = 0.8 # cfl number for hyperbolic system castro.init_shrink = 1.0 # scale back initial timestep castro.change_max = 1.1 # scale back initial timestep -castro.fixed_dt = 5.e-4 +castro.fixed_dt = 1.e-4 # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass @@ -63,7 +63,7 @@ amr.derive_plot_vars = ALL # problem initialization -problem.T0 = 4.e9 +problem.T0 = 5.e9 problem.dT_fact = 0.2 -problem.rho0 = 5.e8 +problem.rho0 = 1.e9 problem.L_pert = 2.e7 diff --git a/Exec/reacting_tests/nse_test/inputs.256 b/Exec/reacting_tests/nse_test/inputs.256 index 28dbf8e3cd..7265950b9e 100644 --- a/Exec/reacting_tests/nse_test/inputs.256 +++ b/Exec/reacting_tests/nse_test/inputs.256 @@ -1,12 +1,12 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- max_step = 15000 -stop_time = 0.08 +stop_time = 0.025 # PROBLEM SIZE & GEOMETRY geometry.is_periodic = 1 1 1 geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical geometry.prob_lo = 0 0 0 -geometry.prob_hi = 1.e8 1.e8 1.e8 +geometry.prob_hi = 2.e7 2.e7 2.e7 amr.n_cell = 256 256 256 @@ -35,7 +35,7 @@ castro.small_temp = 1.e7 castro.cfl = 0.8 # cfl number for hyperbolic system castro.init_shrink = 1.0 # scale back initial timestep castro.change_max = 1.1 # scale back initial timestep -castro.fixed_dt = 2.5e-4 +castro.fixed_dt = 5.e-5 # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass @@ -63,7 +63,7 @@ amr.derive_plot_vars = ALL # problem initialization -problem.T0 = 4.e9 +problem.T0 = 5.e9 problem.dT_fact = 0.2 -problem.rho0 = 5.e8 +problem.rho0 = 1.e9 problem.L_pert = 2.e7 diff --git a/Exec/reacting_tests/nse_test/inputs.256.test b/Exec/reacting_tests/nse_test/inputs.256.test new file mode 100644 index 0000000000..cba5f26d6e --- /dev/null +++ b/Exec/reacting_tests/nse_test/inputs.256.test @@ -0,0 +1,68 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 15000 +stop_time = 0.01 + +# PROBLEM SIZE & GEOMETRY +geometry.is_periodic = 1 1 1 +geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical +geometry.prob_lo = 0 0 0 +geometry.prob_hi = 1.e7 1.e7 1.e7 +amr.n_cell = 256 256 256 + + +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +# 0 = Interior 3 = Symmetry +# 1 = Inflow 4 = SlipWall +# 2 = Outflow 5 = NoSlipWall +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +castro.lo_bc = 0 0 0 +castro.hi_bc = 0 0 0 + +# WHICH PHYSICS +castro.do_hydro = 1 +castro.do_react = 1 + +castro.ppm_type = 1 +castro.ppm_temp_fix = 0 + +castro.use_flattening = 1 + +castro.riemann_solver = 0 + +castro.small_temp = 1.e7 + +# TIME STEP CONTROL +castro.cfl = 0.8 # cfl number for hyperbolic system +castro.init_shrink = 1.0 # scale back initial timestep +castro.change_max = 1.1 # scale back initial timestep + +# DIAGNOSTICS & VERBOSITY +castro.sum_interval = 1 # timesteps between computing mass +castro.v = 1 # verbosity in Castro.cpp +amr.v = 1 # verbosity in Amr.cpp +#amr.grid_log = grdlog # name of grid logging file + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed +amr.ref_ratio = 2 2 2 2 # refinement ratio +amr.regrid_int = 2 2 2 2 # how often to regrid +amr.blocking_factor = 8 # block factor in grid generation +amr.max_grid_size = 64 +amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est + +# CHECKPOINT FILES +amr.checkpoint_files_output = 0 +amr.check_file = nse_test_256_chk # root name of checkpoint file +amr.check_int = 300 # number of timesteps between checkpoints + +# PLOTFILES +amr.plot_file = nse_test_256_plt # root name of plotfile +amr.plot_per = 0.24 +amr.derive_plot_vars = ALL + +# problem initialization + +problem.T0 = 5.e9 +problem.dT_fact = 0.2 +problem.rho0 = 5.e8 +problem.L_pert = 2.e7 diff --git a/Exec/reacting_tests/nse_test/inputs.512 b/Exec/reacting_tests/nse_test/inputs.512 index 4598d1e124..14e56032c3 100644 --- a/Exec/reacting_tests/nse_test/inputs.512 +++ b/Exec/reacting_tests/nse_test/inputs.512 @@ -1,12 +1,12 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- max_step = 15000 -stop_time = 0.08 +stop_time = 0.025 # PROBLEM SIZE & GEOMETRY geometry.is_periodic = 1 1 1 geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical geometry.prob_lo = 0 0 0 -geometry.prob_hi = 1.e8 1.e8 1.e8 +geometry.prob_hi = 2.e7 2.e7 2.e7 amr.n_cell = 512 512 512 @@ -35,7 +35,7 @@ castro.small_temp = 1.e7 castro.cfl = 0.8 # cfl number for hyperbolic system castro.init_shrink = 1.0 # scale back initial timestep castro.change_max = 1.1 # scale back initial timestep -castro.fixed_dt = 1.25e-4 +castro.fixed_dt = 2.5e-5 # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass @@ -63,7 +63,7 @@ amr.derive_plot_vars = ALL # problem initialization -problem.T0 = 4.e9 +problem.T0 = 5.e9 problem.dT_fact = 0.2 -problem.rho0 = 5.e8 +problem.rho0 = 1.e9 problem.L_pert = 2.e7 diff --git a/Exec/reacting_tests/nse_test/inputs.64 b/Exec/reacting_tests/nse_test/inputs.64 index 3bed4d0435..99233c50f3 100644 --- a/Exec/reacting_tests/nse_test/inputs.64 +++ b/Exec/reacting_tests/nse_test/inputs.64 @@ -1,12 +1,12 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- max_step = 15000 -stop_time = 0.08 +stop_time = 0.025 # PROBLEM SIZE & GEOMETRY geometry.is_periodic = 1 1 1 geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical geometry.prob_lo = 0 0 0 -geometry.prob_hi = 1.e8 1.e8 1.e8 +geometry.prob_hi = 2.e7 2.e7 2.e7 amr.n_cell = 64 64 64 @@ -35,7 +35,7 @@ castro.small_temp = 1.e7 castro.cfl = 0.8 # cfl number for hyperbolic system castro.init_shrink = 1.0 # scale back initial timestep castro.change_max = 1.1 # scale back initial timestep -castro.fixed_dt = 1.e-3 +castro.fixed_dt = 2.e-4 # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass @@ -63,7 +63,7 @@ amr.derive_plot_vars = ALL # problem initialization -problem.T0 = 4.e9 +problem.T0 = 5.e9 problem.dT_fact = 0.2 -problem.rho0 = 5.e8 +problem.rho0 = 1.e9 problem.L_pert = 2.e7 diff --git a/Exec/reacting_tests/nse_test/problem_initialize.H b/Exec/reacting_tests/nse_test/problem_initialize.H index 8fd1773a3e..cb885ded7c 100644 --- a/Exec/reacting_tests/nse_test/problem_initialize.H +++ b/Exec/reacting_tests/nse_test/problem_initialize.H @@ -24,8 +24,8 @@ void problem_initialize () // we only work in NSE mode, so the mass fractions are not used by // the EOS -#ifndef NSE_TABLE - amrex::Error("Error: this problem requires USE_NSE=TRUE"); +#if !defined(NSE_TABLE) && !defined(NSE_NET) + amrex::Error("Error: this problem requires USE_NSE_NET=TRUE or USE_NSE_TABLE=TRUE"); #endif } diff --git a/Exec/reacting_tests/nse_test/problem_initialize_state_data.H b/Exec/reacting_tests/nse_test/problem_initialize_state_data.H index f9058fa595..f2956f812b 100644 --- a/Exec/reacting_tests/nse_test/problem_initialize_state_data.H +++ b/Exec/reacting_tests/nse_test/problem_initialize_state_data.H @@ -4,6 +4,11 @@ #include #include #include +#if defined(NSE_TABLE) +#include +#elif defined(NSE_NET) +#include +#endif AMREX_GPU_HOST_DEVICE AMREX_INLINE void problem_initialize_state_data (int i, int j, int k, @@ -49,56 +54,81 @@ void problem_initialize_state_data (int i, int j, int k, ye = ye0; } - state(i,j,k,UMX) = 0.0_rt; - state(i,j,k,UMY) = 0.0_rt; - state(i,j,k,UMZ) = 0.0_rt; + state(i,j,k,UMX) = problem::rho0 * problem::u0; + state(i,j,k,UMY) = problem::rho0 * problem::v0; + state(i,j,k,UMZ) = problem::rho0 * problem::w0; // we are isentropic, so find rho - eos_t eos_state; - eos_state.T = T; - eos_state.rho = problem::rho0; + burn_t burn_state; + burn_state.T = T; + burn_state.rho = problem::rho0; + burn_state.y_e = ye; - Real xn[NumSpec]; +#if defined(NSE_TABLE) Real abar; Real dq; Real dyedt; - nse_interp(T, problem::rho0, ye, abar, dq, dyedt, xn); + nse_interp(T, problem::rho0, ye, abar, dq, dyedt, burn_state.xn); +#elif defined(NSE_NET) + Real eps = 1.e-10_rt; + bool input_ye_is_valid = true; + + // we need an initial guess for the mass fractions -- make them all the same + for (int n = 0; n < NumSpec; ++n) { + burn_state.xn[n] = 1.0_rt / static_cast(NumSpec); + } + + // we also need a guess for the initial chemical potential + burn_state.mu_p = -3.0_rt; + burn_state.mu_n = -12.0_rt; + + auto nse_state = get_actual_nse_state(burn_state, eps, input_ye_is_valid); + + for (int n = 0; n < NumSpec; ++n) { + burn_state.xn[n] = nse_state.xn[n]; + } +#endif // since the species are interpolated, normalize them Real sumX = 0.0_rt; for (int n = 0; n < NumSpec; n++) { - sumX += xn[n]; + sumX += burn_state.xn[n]; } for (int n = 0; n < NumSpec; n++) { - xn[n] /= sumX; + burn_state.xn[n] /= sumX; } - for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = xn[n]; - } - eos_state.aux[AuxZero::iye] = ye; - eos_state.aux[AuxZero::iabar] = abar; - eos_state.aux[AuxZero::ibea] = dq; +#ifdef AUX_THERMO + burn_state.aux[AuxZero::iye] = ye; + burn_state.aux[AuxZero::iabar] = abar; + burn_state.aux[AuxZero::ibea] = dq; +#endif + + eos(eos_input_rt, burn_state); - eos(eos_input_rt, eos_state); + state(i,j,k,URHO) = burn_state.rho; - state(i,j,k,URHO) = eos_state.rho; + state(i,j,k,UEDEN) = burn_state.rho * burn_state.e + + 0.5_rt * burn_state.rho * (problem::u0 * problem::u0 + + problem::v0 * problem::v0 + + problem::w0 * problem::w0); - state(i,j,k,UEDEN) = eos_state.rho * eos_state.e; - state(i,j,k,UEINT) = eos_state.rho * eos_state.e; + state(i,j,k,UEINT) = burn_state.rho * burn_state.e; for (int n = 0; n < NumSpec; n++) { - state(i,j,k,UFS+n) = state(i,j,k,URHO) * xn[n]; + state(i,j,k,UFS+n) = state(i,j,k,URHO) * burn_state.xn[n]; } +#ifdef AUX_THERMO for (int n = 0; n < NumAux; n++) { - state(i,j,k,UFX+n) = state(i,j,k,URHO) * eos_state.aux[n]; + state(i,j,k,UFX+n) = state(i,j,k,URHO) * burn_state.aux[n]; } +#endif - state(i,j,k,UTEMP) = eos_state.T; + state(i,j,k,UTEMP) = burn_state.T; } diff --git a/Exec/reacting_tests/reacting_bubble/Problem_Derive.cpp b/Exec/reacting_tests/reacting_bubble/Problem_Derive.cpp index a3ecdc4019..96ce08233a 100644 --- a/Exec/reacting_tests/reacting_bubble/Problem_Derive.cpp +++ b/Exec/reacting_tests/reacting_bubble/Problem_Derive.cpp @@ -53,7 +53,7 @@ void ca_derpi(const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, for (int n = 0; n < NumSpec; n++) { eos_state.xn[n] = dat(i,j,k,UFS+n) / dat(i,j,k,URHO); } -#if NAUX > 0 +#if NAUX_NET > 0 for (int n = 0; n < NumAux; n++) { eos_state.aux[n] = dat(i,j,k,UFX+n) / dat(i,j,k,URHO); } @@ -112,7 +112,7 @@ void ca_derpioverp0(const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, for (int n = 0; n < NumSpec; n++) { eos_state.xn[n] = dat(i,j,k,UFS+n) / dat(i,j,k,URHO); } -#if NAUX > 0 +#if NAUX_NET > 0 for (int n = 0; n < NumAux; n++) { eos_state.aux[n] = dat(i,j,k,UFX+n) / dat(i,j,k,URHO); } diff --git a/Exec/science/Detonation/GNUmakefile.nse b/Exec/science/Detonation/GNUmakefile.nse index 31275fffd1..c105f3f15d 100644 --- a/Exec/science/Detonation/GNUmakefile.nse +++ b/Exec/science/Detonation/GNUmakefile.nse @@ -19,7 +19,7 @@ EOS_DIR := helmholtz # This sets the EOS directory in $(MICROPHYSICS_HOME)/networks NETWORK_DIR := aprox19 -USE_NSE = TRUE +USE_NSE_TABLE = TRUE Bpack := ./Make.package Blocs := . diff --git a/Exec/science/Detonation/problem_initialize.H b/Exec/science/Detonation/problem_initialize.H index cdb43c88e5..ab60052bd4 100644 --- a/Exec/science/Detonation/problem_initialize.H +++ b/Exec/science/Detonation/problem_initialize.H @@ -59,6 +59,11 @@ void problem_initialize () eos_state.xn[n] = problem::ambient_comp[n]; } +#ifdef AUX_THERMO + // set the aux quantities -- we need to do this if we are using the NSE network + set_aux_comp_from_X(eos_state); +#endif + eos_state.T = problem::T_l; eos(eos_input_rt, eos_state); diff --git a/Exec/science/flame/analysis/flame_speed.py b/Exec/science/flame/analysis/flame_speed.py index e718022272..e91fe3dfc9 100755 --- a/Exec/science/flame/analysis/flame_speed.py +++ b/Exec/science/flame/analysis/flame_speed.py @@ -1,6 +1,5 @@ #!/usr/bin/env python3 -from __future__ import print_function import argparse import matplotlib @@ -17,7 +16,7 @@ class Profile: """read a plotfile using yt and store the 1d profile for T and enuc""" def __init__(self, plotfile): - ds = yt.load(plotfile) + ds = yt.load(plotfile, hint="castro") time = float(ds.current_time) ad = ds.all_data() @@ -108,7 +107,7 @@ def find_flame_width(self): for n in range(0, len(plot_nums), args.skip): - pfile = "{}{}".format(prefix, plot_nums[n]) + pfile = f"{prefix}{plot_nums[n]}" p = Profile(pfile) x1 = p.find_x_for_T(T_0=1.5e9) @@ -158,4 +157,4 @@ def find_flame_width(self): plt.savefig("speed.png") for to, vo, vso, wo in zip(t, v, vs, w): - print("{:10.3g} : {:15.8g} +/- {:15.8g} | {:15.8g}".format(to, vo, vso, wo)) + print(f"{to:10.3g} : {vo:15.8g} +/- {vso:15.8g} | {wo:15.8g}") diff --git a/Exec/science/flame/analysis/flame_speed_convergence.py b/Exec/science/flame/analysis/flame_speed_convergence.py index 3b402c805b..0d42188cf8 100755 --- a/Exec/science/flame/analysis/flame_speed_convergence.py +++ b/Exec/science/flame/analysis/flame_speed_convergence.py @@ -55,7 +55,7 @@ def __str__(self): files = glob.glob("*plt?????") + glob.glob("*plt??????") prefix = files[0].split("plt")[0] + "plt" plot_nums = sorted([p.split("plt")[1] for p in files], key=int) - sorted_files = ["{}{}".format(prefix, q) for q in plot_nums] + sorted_files = [f"{prefix}{q}" for q in plot_nums] # we'll operate on the last 2 files. Compute the width and # flame_speed from these diff --git a/Exec/science/flame/analysis/profile_start_finish.py b/Exec/science/flame/analysis/profile_start_finish.py index f88d4813b0..bc702f43be 100755 --- a/Exec/science/flame/analysis/profile_start_finish.py +++ b/Exec/science/flame/analysis/profile_start_finish.py @@ -1,12 +1,9 @@ #!/usr/bin/env python3 -from __future__ import print_function import argparse -import sys import numpy as np -from cycler import cycler import matplotlib matplotlib.use('agg') import matplotlib.pyplot as plt @@ -15,7 +12,7 @@ def get_Te_profile(plotfile): - ds = yt.load(plotfile) + ds = yt.load(plotfile, hint="castro") time = float(ds.current_time) ad = ds.all_data() @@ -45,14 +42,14 @@ def doit(prefix, nums, xmax): for n in range(0, len(nums)): - pfile = "{}{}".format(prefix, nums[n]) + pfile = f"{prefix}{nums[n]}" time, x, T, p, rho, v, enuc = get_Te_profile(pfile) - ax_T.plot(x, T, color="C{}".format(n), label="t = {:6.4g} s".format(time)) - ax_r.plot(x, rho, color="C{}".format(n)) - ax_p.plot(x, p, color="C{}".format(n)) - ax_v.plot(x, v, color="C{}".format(n)) + ax_T.plot(x, T, color=f"C{n}", label=f"t = {time:6.4g} s") + ax_r.plot(x, rho, color=f"C{n}") + ax_p.plot(x, p, color=f"C{n}") + ax_v.plot(x, v, color=f"C{n}") ax_T.legend(frameon=False, fontsize="small") ax_T.set_ylabel("T (K)") diff --git a/Exec/science/flame/analysis/profiles.py b/Exec/science/flame/analysis/profiles.py index 302f3b285d..db11808090 100755 --- a/Exec/science/flame/analysis/profiles.py +++ b/Exec/science/flame/analysis/profiles.py @@ -22,7 +22,7 @@ def rgba_to_hex(rgba): def get_Te_profile(plotfile): - ds = yt.load(plotfile) + ds = yt.load(plotfile, hint="castro") time = float(ds.current_time) ad = ds.all_data() diff --git a/Exec/science/flame/analysis/profiles_zoom.py b/Exec/science/flame/analysis/profiles_zoom.py index 77dabc6b11..536ac5bf5b 100755 --- a/Exec/science/flame/analysis/profiles_zoom.py +++ b/Exec/science/flame/analysis/profiles_zoom.py @@ -1,6 +1,5 @@ #!/usr/bin/env python3 -from __future__ import print_function import argparse import matplotlib @@ -13,7 +12,7 @@ def get_Te_profile(plotfile): - ds = yt.load(plotfile) + ds = yt.load(plotfile, hint="castro") time = float(ds.current_time) ad = ds.all_data() @@ -35,16 +34,16 @@ def doit(pprefix, nums, skip): ax_T = f.add_subplot(211) ax_e = f.add_subplot(212) - + for n in range(0, len(nums), skip): - pfile = "{}{}".format(prefix, nums[n]) + pfile = f"{prefix}{nums[n]}" time, x, T, enuc = get_Te_profile(pfile) - ax_T.plot(x, T, 'rx', label="t = {:6.4g} s".format(time)) + ax_T.plot(x, T, 'rx', label=f"t = {time:6.4g} s") ax_e.plot(x, enuc) - + ax_T.legend(frameon=False) ax_T.set_ylabel("T (K)") ax_T.set_xlim(22000, 28000) @@ -56,7 +55,7 @@ def doit(pprefix, nums, skip): ax_e.set_xlim(22000, 28000) f.savefig("flame.png") - + if __name__ == "__main__": diff --git a/Exec/science/flame/flame_wave_tests/aprox13/inputs.1d.flame_wave.boost b/Exec/science/flame/flame_wave_tests/aprox13/inputs.1d.flame_wave.boost index 0f39f69e40..f27d1ea6d3 100644 --- a/Exec/science/flame/flame_wave_tests/aprox13/inputs.1d.flame_wave.boost +++ b/Exec/science/flame/flame_wave_tests/aprox13/inputs.1d.flame_wave.boost @@ -40,7 +40,6 @@ castro.use_retry = 1 # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass -amr.data_log = "toy_flame.log" castro.v = 1 # verbosity in Castro.cpp amr.v = 1 # verbosity in Amr.cpp diff --git a/Exec/science/flame/flame_wave_tests/triple_alpha_plus_cago/inputs.1d.flame_wave.boost b/Exec/science/flame/flame_wave_tests/triple_alpha_plus_cago/inputs.1d.flame_wave.boost index 4172f6d4a2..c82289231a 100644 --- a/Exec/science/flame/flame_wave_tests/triple_alpha_plus_cago/inputs.1d.flame_wave.boost +++ b/Exec/science/flame/flame_wave_tests/triple_alpha_plus_cago/inputs.1d.flame_wave.boost @@ -41,7 +41,6 @@ castro.use_retry = 1 # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass -amr.data_log = "toy_flame.log" castro.v = 1 # verbosity in Castro.cpp amr.v = 1 # verbosity in Amr.cpp diff --git a/Exec/science/flame/inputs.1d b/Exec/science/flame/inputs.1d index b6718f837e..0aa281afee 100644 --- a/Exec/science/flame/inputs.1d +++ b/Exec/science/flame/inputs.1d @@ -39,7 +39,6 @@ castro.change_max = 1.1 # max time step growth # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass -amr.data_log = "toy_flame.log" castro.v = 1 # verbosity in Castro.cpp amr.v = 1 # verbosity in Amr.cpp diff --git a/Exec/science/flame/inputs.1d.aprox13 b/Exec/science/flame/inputs.1d.aprox13 index 29b6476026..72fb159b54 100644 --- a/Exec/science/flame/inputs.1d.aprox13 +++ b/Exec/science/flame/inputs.1d.aprox13 @@ -39,7 +39,6 @@ castro.change_max = 1.1 # max time step growth # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass -amr.data_log = "toy_flame.log" castro.v = 1 # verbosity in Castro.cpp amr.v = 1 # verbosity in Amr.cpp diff --git a/Exec/science/flame/inputs.1d.sdc b/Exec/science/flame/inputs.1d.sdc index 2da434f0f0..e20436e736 100644 --- a/Exec/science/flame/inputs.1d.sdc +++ b/Exec/science/flame/inputs.1d.sdc @@ -55,7 +55,6 @@ castro.dtnuc_e = 0.25 # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass -amr.data_log = "toy_flame.log" castro.v = 1 # verbosity in Castro.cpp amr.v = 1 # verbosity in Amr.cpp diff --git a/Exec/science/flame/inputs.1d.sdc.test b/Exec/science/flame/inputs.1d.sdc.test index 7821274d7e..6952901fd3 100644 --- a/Exec/science/flame/inputs.1d.sdc.test +++ b/Exec/science/flame/inputs.1d.sdc.test @@ -42,7 +42,6 @@ castro.change_max = 1.1 # max time step growth # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass -amr.data_log = "toy_flame.log" castro.v = 1 # verbosity in Castro.cpp amr.v = 1 # verbosity in Amr.cpp diff --git a/Exec/science/flame/inputs.1d.sdc.test.3alpha b/Exec/science/flame/inputs.1d.sdc.test.3alpha index d103e5b11d..5cb9b6d23a 100644 --- a/Exec/science/flame/inputs.1d.sdc.test.3alpha +++ b/Exec/science/flame/inputs.1d.sdc.test.3alpha @@ -42,7 +42,6 @@ castro.change_max = 1.1 # max time step growth # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass -amr.data_log = "toy_flame.log" castro.v = 1 # verbosity in Castro.cpp amr.v = 1 # verbosity in Amr.cpp diff --git a/Exec/science/flame/inputs.H_He b/Exec/science/flame/inputs.H_He index ab2d768c6b..65c47fdaf4 100644 --- a/Exec/science/flame/inputs.H_He +++ b/Exec/science/flame/inputs.H_He @@ -39,7 +39,6 @@ castro.change_max = 1.1 # max time step growth # DIAGNOSTICS & VERBOSITY castro.sum_interval = 1 # timesteps between computing mass -amr.data_log = "toy_flame.log" castro.v = 1 # verbosity in Castro.cpp amr.v = 1 # verbosity in Amr.cpp diff --git a/Exec/science/flame/problem_diagnostics.H b/Exec/science/flame/problem_diagnostics.H new file mode 100644 index 0000000000..6408f7a238 --- /dev/null +++ b/Exec/science/flame/problem_diagnostics.H @@ -0,0 +1,70 @@ +#ifndef problem_diagnostics_H +#define problem_diagnostics_H + +AMREX_INLINE +void +Castro::problem_diagnostics () +{ + int finest_level = parent->finestLevel(); + Real time = state[State_Type].curTime(); + + Real T_max = -1.0e200; + Real T_min = 1.0e200; + Real grad_T_max = -1.0e200; + Real rho_fuel_dot = 0.0; + Real rho_fuel_initial = 1.0; + + Real flame_width = 0.0; + Real flame_speed = 0.0; + + for (int lev = 0; lev <= finest_level; lev++) + { + Castro& ca_lev = getLevel(lev); + + ca_lev.flame_width_properties(time, T_max, T_min, grad_T_max); + + // For now we only do the flame speed calculation on the coarse grid since + // it doesn't share the same timestep as the fine grids. + + if (lev == 0) + ca_lev.flame_speed_properties(time, rho_fuel_dot); + } + + ParallelDescriptor::ReduceRealMax(T_max); + ParallelDescriptor::ReduceRealMin(T_min); + ParallelDescriptor::ReduceRealMax(grad_T_max); + + flame_width = (T_max - T_min) / grad_T_max; + + ParallelDescriptor::ReduceRealSum(rho_fuel_dot); + + // Note that rho_fuel_dot has already been multiplied by dx, so + // the dimensionality here checks out. + + flame_speed = -rho_fuel_dot / rho_fuel_initial; + + if (ParallelDescriptor::IOProcessor()) { + + if (verbose > 0) { + std::cout << '\n'; + std::cout << "TIME= " << time << " FLAME WIDTH = " << flame_width << '\n'; + std::cout << "TIME= " << time << " FLAME SPEED = " << flame_speed << '\n'; + } + + std::ostream& datalog = *Castro::problem_data_logs[0]; + + if (time == 0.0) { + datalog << std::setw(14) << " time "; + datalog << std::setw(14) << " flame width "; + datalog << std::setw(14) << " flame speed " << std::endl; + } + + // Write the quantities at this time + datalog << std::setw(14) << time; + datalog << std::setw(14) << std::setprecision(6) << flame_width; + datalog << std::setw(14) << std::setprecision(6) << flame_speed << std::endl; + + } +} + +#endif diff --git a/Exec/science/flame/problem_initialize.H b/Exec/science/flame/problem_initialize.H index 9cf6969137..da36f7dcdb 100644 --- a/Exec/science/flame/problem_initialize.H +++ b/Exec/science/flame/problem_initialize.H @@ -107,5 +107,20 @@ void problem_initialize () // mass flux will be constant across the flame problem::mass_flux = problem::rho_fuel * problem::v_inflow; + + // Set up Castro data logs for this problem + + if (castro::sum_interval > 0 && amrex::ParallelDescriptor::IOProcessor()) { + + Castro::problem_data_logs.resize(1); + + Castro::problem_data_logs[0].reset(new std::fstream); + Castro::problem_data_logs[0]->open("toy_flame.log", std::ios::out | std::ios::app); + if (!Castro::problem_data_logs[0]->good()) { + amrex::FileOpenFailed("toy_flame.log"); + } + + } + } #endif diff --git a/Exec/science/flame/sum_integrated_quantities.cpp b/Exec/science/flame/sum_integrated_quantities.cpp deleted file mode 100644 index e47f824f02..0000000000 --- a/Exec/science/flame/sum_integrated_quantities.cpp +++ /dev/null @@ -1,164 +0,0 @@ -#include - -#include -#include - -using namespace amrex; - -void -Castro::sum_integrated_quantities () -{ - if (verbose <= 0) return; - - bool local_flag = true; - - int finest_level = parent->finestLevel(); - Real time = state[State_Type].curTime(); - Real prev_time = state[State_Type].prevTime(); - Real mass = 0.0; - Real mom[3] = { 0.0 }; - Real com[3] = { 0.0 }; - Real com_vel[3] = { 0.0 }; - Real rho_e = 0.0; - Real rho_K = 0.0; - Real rho_E = 0.0; - Real dt_crse = parent->dtLevel(0); - - Real T_max = -1.0e200; - Real T_min = 1.0e200; - Real grad_T_max = -1.0e200; - Real rho_fuel_dot = 0.0; - Real rho_fuel_initial = 1.0; - - Real flame_width = 0.0; - Real flame_speed = 0.0; - - for (int lev = 0; lev <= finest_level; lev++) - { - Castro& ca_lev = getLevel(lev); - - mass += ca_lev.volWgtSum("density", time, local_flag); - mom[0] += ca_lev.volWgtSum("xmom", time, local_flag); - mom[1] += ca_lev.volWgtSum("ymom", time, local_flag); - mom[2] += ca_lev.volWgtSum("zmom", time, local_flag); - - com[0] += ca_lev.locWgtSum("density", time, 0, local_flag); - com[1] += ca_lev.locWgtSum("density", time, 1, local_flag); - com[2] += ca_lev.locWgtSum("density", time, 2, local_flag); - - rho_e += ca_lev.volWgtSum("rho_e", time, local_flag); - rho_K += ca_lev.volWgtSum("kineng", time, local_flag); - rho_E += ca_lev.volWgtSum("rho_E", time, local_flag); - - - ca_lev.flame_width_properties(time, T_max, T_min, grad_T_max); - - // For now we only do the flame speed calculation on the coarse grid since - // it doesn't share the same timestep as the fine grids. - - if (lev == 0) - ca_lev.flame_speed_properties(time, rho_fuel_dot); - } - - if (verbose > 0) - { - const int nfoo = 7; - Real foo[nfoo] = {mass, mom[0], mom[1], mom[2], rho_e, rho_K, rho_E}; - -#ifdef BL_LAZY - Lazy::QueueReduction( [=] () mutable { -#endif - - ParallelDescriptor::ReduceRealSum(foo, nfoo, ParallelDescriptor::IOProcessorNumber()); - - ParallelDescriptor::ReduceRealSum(com, 3, ParallelDescriptor::IOProcessorNumber()); - - ParallelDescriptor::ReduceRealMax(T_max); - ParallelDescriptor::ReduceRealMin(T_min); - ParallelDescriptor::ReduceRealMax(grad_T_max); - - flame_width = (T_max - T_min) / grad_T_max; - - ParallelDescriptor::ReduceRealSum(rho_fuel_dot); - - // Note that rho_fuel_dot has already been multiplied by dx, so - // the dimensionality here checks out. - - flame_speed = -rho_fuel_dot / rho_fuel_initial; - - if (ParallelDescriptor::IOProcessor()) { - - int i = 0; - mass = foo[i++]; - mom[0] = foo[i++]; - mom[1] = foo[i++]; - mom[2] = foo[i++]; - rho_e = foo[i++]; - rho_K = foo[i++]; - rho_E = foo[i++]; - - std::cout << '\n'; - std::cout << "TIME= " << time << " MASS = " << mass << '\n'; - std::cout << "TIME= " << time << " XMOM = " << mom[0] << '\n'; - std::cout << "TIME= " << time << " YMOM = " << mom[1] << '\n'; - std::cout << "TIME= " << time << " ZMOM = " << mom[2] << '\n'; - std::cout << "TIME= " << time << " RHO*e = " << rho_e << '\n'; - std::cout << "TIME= " << time << " RHO*K = " << rho_K << '\n'; - std::cout << "TIME= " << time << " RHO*E = " << rho_E << '\n'; - std::cout << "TIME= " << time << " FLAME WIDTH = " << flame_width << '\n'; - std::cout << "TIME= " << time << " FLAME SPEED = " << flame_speed << '\n'; - - if (parent->NumDataLogs() > 0 ) { - - std::ostream& data_log1 = parent->DataLog(0); - - if (data_log1.good()) { - - if (time == 0.0) { - data_log1 << std::setw(14) << " time "; - data_log1 << std::setw(14) << " mass "; - data_log1 << std::setw(14) << " xmom "; - data_log1 << std::setw(14) << " ymom "; - data_log1 << std::setw(14) << " zmom "; - data_log1 << std::setw(14) << " rho_K "; - data_log1 << std::setw(14) << " rho_e "; - data_log1 << std::setw(14) << " rho_E "; - data_log1 << std::setw(14) << " flame width "; - data_log1 << std::setw(14) << " flame speed " << std::endl; - } - - // Write the quantities at this time - data_log1 << std::setw(14) << time; - data_log1 << std::setw(14) << std::setprecision(6) << mass; - data_log1 << std::setw(14) << std::setprecision(6) << mom[0]; - data_log1 << std::setw(14) << std::setprecision(6) << mom[1]; - data_log1 << std::setw(14) << std::setprecision(6) << mom[2]; - data_log1 << std::setw(14) << std::setprecision(6) << rho_K; - data_log1 << std::setw(14) << std::setprecision(6) << rho_e; - data_log1 << std::setw(14) << std::setprecision(6) << rho_E; - data_log1 << std::setw(14) << std::setprecision(6) << flame_width; - data_log1 << std::setw(14) << std::setprecision(6) << flame_speed << std::endl; - - } - - } - - for (int i = 0; i <= 2; i++) { - com[i] = com[i] / mass; - com_vel[i] = mom[i] / mass; - } - - std::cout << "TIME= " << time << " CENTER OF MASS X-LOC = " << com[0] << '\n'; - std::cout << "TIME= " << time << " CENTER OF MASS X-VEL = " << com_vel[0] << '\n'; - - std::cout << "TIME= " << time << " CENTER OF MASS Y-LOC = " << com[1] << '\n'; - std::cout << "TIME= " << time << " CENTER OF MASS Y-VEL = " << com_vel[1] << '\n'; - - std::cout << "TIME= " << time << " CENTER OF MASS Z-LOC = " << com[2] << '\n'; - std::cout << "TIME= " << time << " CENTER OF MASS Z-VEL = " << com_vel[2] << '\n'; - } -#ifdef BL_LAZY - }); -#endif - } -} diff --git a/Exec/science/flame_wave/ci-benchmarks/grid_diag.out b/Exec/science/flame_wave/ci-benchmarks/grid_diag.out index 58b8989715..00b810f9dd 100644 --- a/Exec/science/flame_wave/ci-benchmarks/grid_diag.out +++ b/Exec/science/flame_wave/ci-benchmarks/grid_diag.out @@ -1,5 +1,5 @@ # COLUMN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 # TIMESTEP TIME MASS XMOM YMOM ZMOM ANG. MOM. X ANG. MOM. Y ANG. MOM. Z KIN. ENERGY INT. ENERGY GAS ENERGY GRAV. ENERGY TOTAL ENERGY CENTER OF MASS X-LOC CENTER OF MASS Y-LOC CENTER OF MASS Z-LOC CENTER OF MASS X-VEL CENTER OF MASS Y-VEL CENTER OF MASS Z-VEL MAXIMUM TEMPERATURE MAXIMUM DENSITY MAXIMUM T_S / T_E 0 0.0000000000000000 1.7166788475340266e+20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 3.5939094484738719e+37 3.5939094484738719e+37 0.0000000000000000e+00 3.5939094484738719e+37 2.7306358373068455e+04 7.1254231429487743e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.3624032936555479e+09 3.4338407089146405e+07 0.0000000000000000e+00 - 1 5.8556067872049896e-08 1.7166803365487742e+20 3.1513080056688146e+23 9.4393200692827477e+22 1.1594142580106661e+20 3.0429100322504389e+23 -2.4950450852898835e+24 2.1763212263548750e+27 7.5396363812524675e+29 3.5939161397498687e+37 3.5939162151462194e+37 0.0000000000000000e+00 3.5939162151462194e+37 2.7306358430778866e+04 7.1254177204721407e+02 0.0000000000000000e+00 1.8356987836210822e+03 5.4985892645916954e+02 6.7538156832480556e-01 1.3624647489781995e+09 3.4338736190485895e+07 3.6445955142666982e-04 - 2 1.2296774253130479e-07 1.7166821577112974e+20 6.6177191170810717e+23 2.0512769223276993e+23 5.1128643455286516e+20 1.3418692305321693e+24 -1.1002931760322212e+25 4.7588171814049066e+27 2.8702315946888442e+30 3.5939240430623659e+37 3.5939243300855208e+37 0.0000000000000000e+00 3.5939243300855208e+37 2.7306358617303587e+04 7.1254115608754887e+02 0.0000000000000000e+00 3.8549472232552935e+03 1.1949078127907428e+03 2.9783407036426519e+00 1.3628615188759754e+09 3.4339036354305416e+07 3.6432541934204224e-04 + 1 5.8556067872049896e-08 1.7166803365487749e+20 3.1513080056552834e+23 9.4393200694787844e+22 1.1594142580066235e+20 3.0429100322420758e+23 -2.4950450852816581e+24 2.1763212264094320e+27 7.5396363801382249e+29 3.5939161397498158e+37 3.5939162151461684e+37 0.0000000000000000e+00 3.5939162151461684e+37 2.7306358430778862e+04 7.1254177204721316e+02 0.0000000000000000e+00 1.8356987836131991e+03 5.4985892647058881e+02 6.7538156832245033e-01 1.3624647490176616e+09 3.4338736190485746e+07 3.6445955141358745e-04 + 2 1.2296774253130479e-07 1.7166821577112987e+20 6.6177191170323842e+23 2.0512769235398071e+23 5.1128643454995517e+20 1.3418692305260181e+24 -1.1002931760262338e+25 4.7588171847234185e+27 2.8702315944527227e+30 3.5939240430622091e+37 3.5939243300853715e+37 0.0000000000000000e+00 3.5939243300853715e+37 2.7306358617303569e+04 7.1254115608754830e+02 0.0000000000000000e+00 3.8549472232269291e+03 1.1949078134968177e+03 2.9783407036256984e+00 1.3628615188920422e+09 3.4339036354304597e+07 3.6432541942946792e-04 diff --git a/Exec/science/flame_wave/ci-benchmarks/species_diag.out b/Exec/science/flame_wave/ci-benchmarks/species_diag.out index 7bed8b6201..4f7ad6aa71 100644 --- a/Exec/science/flame_wave/ci-benchmarks/species_diag.out +++ b/Exec/science/flame_wave/ci-benchmarks/species_diag.out @@ -1,5 +1,5 @@ # COLUMN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 # TIMESTEP TIME Mass He4 Mass C12 Mass O16 Mass Ne20 Mass Mg24 Mass Si28 Mass S32 Mass Ar36 Mass Ca40 Mass Ti44 Mass Cr48 Mass Fe52 Mass Ni56 0 0.0000000000000000 2.5645624865668260e-15 8.6334683437635802e-24 8.6334683437635802e-24 8.6334683437635802e-24 8.6334683437635802e-24 8.6334683437635802e-24 8.6334683437635802e-24 8.6334683437635802e-24 8.6334683437635802e-24 8.6334683437635802e-24 8.6334683437635802e-24 8.6334683437635802e-24 8.3770120959700128e-14 - 1 5.8556067872049896e-08 2.5645578286502653e-15 4.6665271222083047e-21 8.6558492055718780e-24 8.6336424231163498e-24 8.6337819380890799e-24 8.6336288308871613e-24 8.6334771302909650e-24 8.6334762294408893e-24 8.6334758774055529e-24 8.6334758341207241e-24 8.6334758323582423e-24 8.6334758322823729e-24 8.3770195844770837e-14 - 2 1.2296774253130479e-07 2.5645527047237343e-15 9.7903860219820578e-21 8.7222124044551565e-24 8.6345521905800748e-24 8.6341373318324136e-24 8.6338084562916837e-24 8.6334877298453766e-24 8.6334858256067090e-24 8.6334850860334561e-24 8.6334849950921233e-24 8.6334849913892191e-24 8.6334849912298162e-24 8.3770287434115253e-14 + 1 5.8556067872049896e-08 2.5645578286513752e-15 4.6665255832439506e-21 8.6562660805528712e-24 8.6336501283765672e-24 8.6337820442917443e-24 8.6336288504892702e-24 8.6334771304039403e-24 8.6334762294436341e-24 8.6334758774059056e-24 8.6334758341207285e-24 8.6334758323582423e-24 8.6334758322823700e-24 8.3770195844770875e-14 + 2 1.2296774253130479e-07 2.5645527047252351e-15 9.7903839457371030e-21 8.7227655843562583e-24 8.6345746607885893e-24 8.6341379740187793e-24 8.6338084924592057e-24 8.6334877300430889e-24 8.6334858256111406e-24 8.6334850860340439e-24 8.6334849950921394e-24 8.6334849913892147e-24 8.6334849912298103e-24 8.3770287434115241e-14 diff --git a/Exec/science/massive_star/GNUmakefile b/Exec/science/massive_star/GNUmakefile index df2652a885..638579dc93 100644 --- a/Exec/science/massive_star/GNUmakefile +++ b/Exec/science/massive_star/GNUmakefile @@ -13,7 +13,7 @@ USE_MPI = TRUE USE_GRAV = TRUE USE_REACT = TRUE -USE_NSE = TRUE +USE_NSE_TABLE = TRUE BL_NO_FORT := TRUE diff --git a/Exec/science/massive_star/problem_initialize_state_data.H b/Exec/science/massive_star/problem_initialize_state_data.H index 610e715688..b5a8bfe433 100644 --- a/Exec/science/massive_star/problem_initialize_state_data.H +++ b/Exec/science/massive_star/problem_initialize_state_data.H @@ -5,7 +5,7 @@ #include #include #ifdef NSE_TABLE -#include +#include #endif AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -50,19 +50,32 @@ void problem_initialize_state_data (int i, int j, int k, state(i,j,k,UFX+AuxZero::iye) = interpolate(dist, model::iaux+AuxZero::iye); - eos_t eos_state; - eos_state.rho = state(i,j,k,URHO); - eos_state.T = state(i,j,k,UTEMP); + burn_t burn_state; + burn_state.rho = state(i,j,k,URHO); + burn_state.T = state(i,j,k,UTEMP); for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = state(i,j,k,UFS+n); + burn_state.xn[n] = state(i,j,k,UFS+n); } - eos_state.aux[AuxZero::iye] = state(i,j,k,UFX+AuxZero::iye); + burn_state.aux[AuxZero::iye] = state(i,j,k,UFX+AuxZero::iye); + +#ifdef SIMPLIFIED_SDC + // if we are doing simplified-SDC + NSE, then the `in_nse()` + // check will use burn_state.y[], so we need to ensure that + // those are initialized + for (int n = 0; n < NumSpec; ++n) { + burn_state.y[SFS+n] = burn_state.rho * burn_state.xn[n]; + } + + // not used + burn_state.y[SEINT] = 0.0; +#endif + // now we use Ye if we are in NSE and use X if we are not. In // either case, we sync up the other composition variables. - auto nse_check = in_nse(eos_state); + auto nse_check = in_nse(burn_state); if (nse_check) { @@ -91,31 +104,33 @@ void problem_initialize_state_data (int i, int j, int k, } else { // we are not in NSE, so set the aux quantities from X-- this - // will fill eos_state.aux[] + // will fill burn_state.aux[] - set_aux_comp_from_X(eos_state); + set_aux_comp_from_X(burn_state); for (int n = 0; n < NumAux; n++) { - state(i,j,k,UFX+n) = eos_state.aux[n]; + state(i,j,k,UFX+n) = burn_state.aux[n]; } } // now call the EOS -- reload the eos_state since composition may // have changed - eos_state.rho = state(i,j,k,URHO); - eos_state.T = state(i,j,k,UTEMP); + burn_state.rho = state(i,j,k,URHO); + burn_state.T = state(i,j,k,UTEMP); for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = state(i,j,k,UFS+n); + burn_state.xn[n] = state(i,j,k,UFS+n); } for (int n = 0; n < NumAux; n++) { - eos_state.aux[n] = state(i,j,k,UFX+n); + burn_state.aux[n] = state(i,j,k,UFX+n); } - eos(eos_input_rt, eos_state); + eos(eos_input_rt, burn_state); + + state(i,j,k,UEINT) = state(i,j,k,URHO) * burn_state.e; + state(i,j,k,UEDEN) = state(i,j,k,URHO) * burn_state.e; - state(i,j,k,UEINT) = state(i,j,k,URHO) * eos_state.e; - state(i,j,k,UEDEN) = state(i,j,k,URHO) * eos_state.e; + // add density scaling to composition for (int n = 0; n < NumSpec; n++) { state(i,j,k,UFS+n) = state(i,j,k,URHO) * state(i,j,k,UFS+n); diff --git a/Exec/science/nova/inputs_nova_t7 b/Exec/science/nova/inputs_nova_t7 index c1e8bb5eb5..3441ad864a 100644 --- a/Exec/science/nova/inputs_nova_t7 +++ b/Exec/science/nova/inputs_nova_t7 @@ -1,12 +1,12 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 100000 +max_step = 200000 # PROBLEM SIZE & GEOMETRY geometry.is_periodic = 1 0 geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical geometry.prob_lo = 0 0 -geometry.prob_hi = 12.192e7 6.096e7 -amr.n_cell = 256 128 +geometry.prob_hi = 12.288e7 6.144e7 +amr.n_cell = 384 192 # >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< # 0 = Interior 3 = Symmetry @@ -63,8 +63,8 @@ amr.v = 1 # verbosity in Amr.cpp amr.max_level = 5 # maximum level number allowed amr.ref_ratio = 2 2 2 2 2 2 # refinement ratio amr.regrid_int = 2 2 2 2 2 2 # how often to regrid -amr.blocking_factor = 128 # block factor in grid generation -amr.max_grid_size = 256 +amr.blocking_factor = 32 # block factor in grid generation +#amr.max_grid_size = 256 amr.n_error_buf = 2 2 2 2 2 2 # number of buffer cells in error est # CHECKPOINT FILES diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 9c7efc7af2..9eb70c227f 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -493,7 +493,7 @@ Castro::read_params () // SCF initial model construction can only be done if both // rotation and gravity have been compiled in. -#if (!defined(GRAVITY) || !defined(ROTATION)) +#if !defined(GRAVITY) || !defined(ROTATION) if (do_scf_initial_model) { amrex::Error("SCF initial model construction is only permitted if USE_GRAV=TRUE and USE_ROTATION=TRUE at compile time."); } diff --git a/Source/driver/Derive.cpp b/Source/driver/Derive.cpp index 0fbfc39016..014853f509 100644 --- a/Source/driver/Derive.cpp +++ b/Source/driver/Derive.cpp @@ -1202,21 +1202,34 @@ extern "C" { Real rhoInv = 1.0_rt / dat(i,j,k,URHO); - eos_t eos_state; - eos_state.rho = dat(i,j,k,URHO); - eos_state.T = dat(i,j,k,UTEMP); - eos_state.e = dat(i,j,k,UEINT) * rhoInv; + burn_t burn_state; + burn_state.rho = dat(i,j,k,URHO); + burn_state.T = dat(i,j,k,UTEMP); + burn_state.e = dat(i,j,k,UEINT) * rhoInv; for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = dat(i,j,k,UFS+n) * rhoInv; + burn_state.xn[n] = dat(i,j,k,UFS+n) * rhoInv; } #if NAUX_NET > 0 for (int n = 0; n < NumAux; n++) { - eos_state.aux[n] = dat(i,j,k,UFX+n) * rhoInv; + burn_state.aux[n] = dat(i,j,k,UFX+n) * rhoInv; } #endif - eos(eos_input_re, eos_state); - der(i,j,k,0) = in_nse(eos_state); + eos(eos_input_re, burn_state); + +#ifdef SIMPLIFIED_SDC + // if we are doing simplified-SDC + NSE, then the `in_nse()` + // check will use burn_state.y[], so we need to ensure that + // those are initialized + for (int n = 0; n < NumSpec; ++n) { + burn_state.y[SFS+n] = burn_state.rho * burn_state.xn[n]; + } + + burn_state.y[SEINT] = burn_state.rho * burn_state.e; +#endif + + der(i,j,k,0) = in_nse(burn_state); + }); } #endif diff --git a/Source/driver/sum_integrated_quantities.cpp b/Source/driver/sum_integrated_quantities.cpp index 5a86e4f28c..7ac3cd7891 100644 --- a/Source/driver/sum_integrated_quantities.cpp +++ b/Source/driver/sum_integrated_quantities.cpp @@ -305,7 +305,7 @@ Castro::sum_integrated_quantities () std::cout << "TIME= " << time << " MAXIMUM TEMPERATURE = " << T_max << '\n'; std::cout << "TIME= " << time << " MAXIMUM DENSITY = " << rho_max << '\n'; -#ifdef REACTION +#ifdef REACTIONS std::cout << "TIME= " << time << " MAXIMUM T_S / T_E = " << ts_te_max << '\n'; #endif diff --git a/Source/driver/timestep.cpp b/Source/driver/timestep.cpp index 9cd65d0123..753e2a3b20 100644 --- a/Source/driver/timestep.cpp +++ b/Source/driver/timestep.cpp @@ -99,7 +99,7 @@ Castro::estdt_cfl (int is_new) } else { // method of lines-style constraint is tougher Real dt_tmp = 1.0_rt/dt1; -#if AMREX_SPACEIM >= 2 +#if AMREX_SPACEDIM >= 2 dt_tmp += 1.0_rt/dt2; #endif #if AMREX_SPACEDIM == 3 @@ -315,6 +315,7 @@ Castro::estdt_burning (int is_new) return {ValLocPair{1.e200_rt, idx}}; } + const auto dx = geom.CellSizeArray(); MultiFab& stateMF = is_new ? get_new_data(State_Type) : get_old_data(State_Type); @@ -360,6 +361,12 @@ Castro::estdt_burning (int is_new) burn_t burn_state; +#if AMREX_SPACEDIM == 1 + burn_state.dx = dx[0]; +#else + burn_state.dx = amrex::min(D_DECL(dx[0], dx[1], dx[2])); +#endif + burn_state.rho = S(i,j,k,URHO); burn_state.T = S(i,j,k,UTEMP); burn_state.e = S(i,j,k,UEINT) * rhoInv; @@ -413,14 +420,20 @@ Castro::estdt_burning (int is_new) Real dt_tmp = 1.e200_rt; #ifdef NSE - // we need to use the eos_state interface here because for - // SDC, if we come in with a burn_t, it expects to - // evaluate the NSE criterion based on the conserved state. - eos_t eos_state; - burn_to_eos(burn_state, eos_state); +#ifdef SIMPLIFIED_SDC + // if we are doing simplified-SDC + NSE, then the `in_nse()` + // check will use burn_state.y[], so we need to ensure that + // those are initialized + for (int n = 0; n < NumSpec; ++n) { + burn_state.y[SFS+n] = burn_state.rho * burn_state.xn[n]; + } + + burn_state.y[SEINT] = burn_state.rho * burn_state.e; + +#endif - if (!in_nse(eos_state)) { + if (!in_nse(burn_state)) { #endif dt_tmp = dtnuc_e * e / dedt; #ifdef NSE diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index 9dd2c4cb53..3953bdd7e4 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -81,6 +81,12 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra burn_t burn_state; +#if AMREX_SPACEDIM == 1 + burn_state.dx = dx[0]; +#else + burn_state.dx = amrex::min(D_DECL(dx[0], dx[1], dx[2])); +#endif + // Initialize some data for later. bool do_burn = true; @@ -98,8 +104,14 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra Real rhoInv = 1.0_rt / U(i,j,k,URHO); burn_state.rho = U(i,j,k,URHO); + + // this T is consistent with UEINT because we did an EOS call before + // calling this function + burn_state.T = U(i,j,k,UTEMP); + burn_state.T_fixed = -1.e30_rt; + #ifdef CXX_MODEL_PARSER if (drive_initial_convection) { Real rr[3] = {0.0_rt}; @@ -120,14 +132,11 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra dist = std::sqrt(rr[0] * rr[0] + rr[1] * rr[1] + rr[2] * rr[2]); } - burn_state.T = interpolate(dist, model::itemp); + burn_state.T_fixed = interpolate(dist, model::itemp); } #endif - - burn_state.e = 0.0_rt; // Energy generated by the burn - for (int n = 0; n < NumSpec; ++n) { burn_state.xn[n] = U(i,j,k,UFS+n) * rhoInv; } @@ -182,7 +191,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra if (reactions.contains(i,j,k)) { - reactions(i,j,k,0) = U(i,j,k,URHO) * burn_state.e / dt; + reactions(i,j,k,0) = (U(i,j,k,URHO) * burn_state.e - U(i,j,k,UEINT)) / dt; if (store_omegadot == 1) { if (reactions.contains(i,j,k)) { @@ -219,8 +228,9 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra U(i,j,k,UFX+n) = U(i,j,k,URHO) * burn_state.aux[n]; } #endif - U(i,j,k,UEINT) += U(i,j,k,URHO) * burn_state.e; - U(i,j,k,UEDEN) += U(i,j,k,URHO) * burn_state.e; + Real reint_old = U(i,j,k,UEINT); + U(i,j,k,UEINT) = U(i,j,k,URHO) * burn_state.e; + U(i,j,k,UEDEN) += U(i,j,k,UEINT) - reint_old; } else { @@ -368,6 +378,12 @@ Castro::react_state(Real time, Real dt) burn_t burn_state; +#if AMREX_SPACEDIM == 1 + burn_state.dx = dx[0]; +#else + burn_state.dx = amrex::min(D_DECL(dx[0], dx[1], dx[2])); +#endif + // Initialize some data for later. bool do_burn = true; diff --git a/Source/reactions/_parameters b/Source/reactions/_parameters new file mode 100644 index 0000000000..e15f939404 --- /dev/null +++ b/Source/reactions/_parameters @@ -0,0 +1,4 @@ +@namespace: integrator + +# Do not subtract the initial energy in the integrator +subtract_internal_energy logical .false. 100 diff --git a/Util/scripts/parse_stdout.py b/Util/scripts/parse_stdout.py new file mode 100644 index 0000000000..597120be60 --- /dev/null +++ b/Util/scripts/parse_stdout.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python + +import re +import sys + +RE_STEP = r"^STEP = (\d+)" +RE_LEVEL = r"^\[Level (\d+) step (\d+)\]" +RE_SDC_ITER = r"^Beginning SDC iteration (\d+) of (\d+)" +RE_VODE = r"DVODE:" + +ofile = sys.argv[-1] + +with open(ofile, "r") as of: + + step = 0 + level = 0 + sdc_iter = 0 + + while (line := of.readline()): + + if g := re.search(RE_STEP, line): + step = g.groups()[0] + continue + + if g := re.search(RE_LEVEL, line): + level = g.groups()[0] + continue + + if g := re.search(RE_SDC_ITER, line): + sdc_iter = g.groups()[0] + continue + + if g := re.search(RE_VODE, line): + print(f"step: {step}, level: {level}, sdc iter: {sdc_iter} ", line) + diff --git a/external/Microphysics b/external/Microphysics index 380d048eea..c750b283ab 160000 --- a/external/Microphysics +++ b/external/Microphysics @@ -1 +1 @@ -Subproject commit 380d048eea6580799f0458cf38fe53ffd86de1e8 +Subproject commit c750b283ab868f319a92bf8d7b4da82dfc89644b diff --git a/external/amrex b/external/amrex index 692059ae91..eefe246388 160000 --- a/external/amrex +++ b/external/amrex @@ -1 +1 @@ -Subproject commit 692059ae91b89245e8b8d28de117df10ddfa8ed6 +Subproject commit eefe2463884081a27025973d4a99dc909657b4f0