diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a667afd0..58278d0a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -84,7 +84,7 @@ jobs: runs-on: ${{matrix.os}} strategy: matrix: - os: [ubuntu-24.04, macos-latest] + os: [ubuntu-22.04, macos-latest] build_type: [Release, Debug] enable_eb: [EB-OFF, EB-ON] include: @@ -93,7 +93,7 @@ jobs: comp: llvm procs: $(sysctl -n hw.ncpu) ccache_cache: /Users/runner/Library/Caches/ccache - - os: ubuntu-24.04 + - os: ubuntu-22.04 install_deps: sudo apt-get update && sudo apt-get install mpich libmpich-dev comp: gnu procs: $(nproc) diff --git a/Exec/RegTests/EB_ODEQty/CMakeLists.txt b/Exec/RegTests/EB_ODEQty/CMakeLists.txt index 06198d3c..5ae0ae30 100644 --- a/Exec/RegTests/EB_ODEQty/CMakeLists.txt +++ b/Exec/RegTests/EB_ODEQty/CMakeLists.txt @@ -1,6 +1,6 @@ set(PELE_PHYSICS_EOS_MODEL Fuego) -set(PELE_PHYSICS_CHEMISTRY_MODEL air) -set(PELE_PHYSICS_TRANSPORT_MODEL Constant) +set(PELE_PHYSICS_CHEMISTRY_MODEL Davis) +set(PELE_PHYSICS_TRANSPORT_MODEL Simple) set(PELE_PHYSICS_ENABLE_SPRAY OFF) set(PELE_PHYSICS_SPRAY_FUEL_NUM 0) set(PELE_PHYSICS_ENABLE_SOOT OFF) diff --git a/Exec/RegTests/EB_ODEQty/GNUmakefile b/Exec/RegTests/EB_ODEQty/GNUmakefile index acda3985..4c90fa29 100644 --- a/Exec/RegTests/EB_ODEQty/GNUmakefile +++ b/Exec/RegTests/EB_ODEQty/GNUmakefile @@ -30,9 +30,9 @@ CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp PELELM_NUM_ODE = 3 # PelePhysics -Chemistry_Model = air +Chemistry_Model = Davis Eos_Model = Fuego -Transport_Model = Constant +Transport_Model = Simple PELE_HOME ?= ../../.. include $(PELE_HOME)/Exec/Make.PeleLMeX diff --git a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp index 9e3dce02..d3422995 100644 --- a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp @@ -25,7 +25,7 @@ set_ode_names(Vector& a_ode_names) void problem_modify_ext_sources( - Real /*time*/, + Real time, Real /*dt*/, const MultiFab& state_old, const MultiFab& /*state_new*/, @@ -45,10 +45,27 @@ problem_modify_ext_sources( ParallelFor( *ext_src, [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { - for (int n = 0; n < NUM_ODE; n++) { - Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n); - Real src = prob_parm.ode_srcstrength * pow(10.0, n + 1) * B_n; - ext_src_arr[box_no](i, j, k, FIRSTODE + n) += src; + if (prob_parm.ode_qty_test) { + for (int n = 0; n < NUM_ODE; n++) { + // Source terms for ODE qty test + + Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n); + Real src = prob_parm.ode_srcstrength * pow(10.0, n + 1) * B_n; + ext_src_arr[box_no](i, j, k, FIRSTODE + n) += src; + } + } + + // Source terms for composition test + if (prob_parm.composition_test) { + Real src = 0.0; + + if (time >= prob_parm.extRhoYCO2_ts) { + Real CO2 = state_old_arr[box_no](i, j, k, FIRSTSPEC + CO2_ID); + src = prob_parm.extRhoYCO2 * CO2; + } + + ext_src_arr[box_no](i, j, k, FIRSTSPEC + CO2_ID) += src; + ext_src_arr[box_no](i, j, k, DENSITY) += src; } }); Gpu::streamSynchronize(); diff --git a/Exec/RegTests/EB_ODEQty/README.md b/Exec/RegTests/EB_ODEQty/README.md index c04720a5..a46ee30e 100644 --- a/Exec/RegTests/EB_ODEQty/README.md +++ b/Exec/RegTests/EB_ODEQty/README.md @@ -1,2 +1,4 @@ ## EB\_ODEQty -A 2D inflow/outflow setup with an optional EB cylinder in the middle of the flow. Demonstrates how to use ProblemSpecificFunctions and the ODE quantities. The ODE quantities experience simple exponential decay that gets stiffer for each quantity. Specifically, $\frac{\partial B_k}{\partial t} = \gamma \cdot 10^{k+1} B_k$, for $k = 0, 1,\dots,$ NUM_ODE and $\gamma < 0$. +This RegTest contains two test cases that demonstrate the use of user defined external source terms: +1. `prob_parm.ode_qty_test = true`: A 2D inflow/outflow case with an optional EB cylinder in the middle of the flow. This demonstrates how to use `ProblemSpecificFunctions` and the ODE quantities. The ODE quantities experience simple exponential decay that gets stiffer for each quantity. Specifically, $\frac{\partial B_k}{\partial t} = \gamma \cdot 10^{k+1} B_k$, for $k = 0, 1,\dots,$ NUM_ODE and $\gamma < 0$. +2. `prob_parm.composition_test = true`: This case tests the implementation of the user defined external sources for density, $\rho$, and the species $Y_m$. The flow is set to zero with initial mass fractions $Y_{CO_2} = 0.15$, $Y_{N_2} = 0.25$ and $Y_{AR} = 0.6$. After $t_s$ seconds, CO2 is added via an external source at a rate of $S_{ext,\rho Y_{CO_2}}$ kg/m$^3$/s. Since CO2 is being added to the domain, density will increase at a rate of $S_{ext,\rho} = S_{ext,\rho Y_{CO_2}}$. Results are compared to analytical solutions in `plot-comps.py` and tested in the CI via `test.py`. diff --git a/Exec/RegTests/EB_ODEQty/composition-test-2d.inp b/Exec/RegTests/EB_ODEQty/composition-test-2d.inp new file mode 100644 index 00000000..0e6b60e6 --- /dev/null +++ b/Exec/RegTests/EB_ODEQty/composition-test-2d.inp @@ -0,0 +1,58 @@ +#----------------------DOMAIN DEFINITION------------------------- +geometry.is_periodic = 1 1 # For each dir, 0: non-perio, 1: periodic +geometry.coord_sys = 0 # 0 => cart, 1 => RZ +geometry.prob_lo = 0 0 +geometry.prob_hi = 0.01 0.01 + +#--------------------------BC FLAGS------------------------------ +# Interior, Inflow, Outflow, Symmetry, +# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm +peleLM.lo_bc = Interior Interior +peleLM.hi_bc = Interior Interior + +#-------------------------AMR CONTROL---------------------------- +amr.n_cell = 1 1 # Level 0 number of cells in each direction +amr.v = 0 # AMR verbose +amr.max_level = 0 # maximum level number allowed +amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est +amr.grid_eff = 0.7 # what constitutes an efficient grid +amr.blocking_factor = 1 # block factor in grid generation (min box size) +amr.max_grid_size = 64 # max box size + +#--------------------------- Problem ---------------------------- +prob.T_mean = 300.0 +prob.P_mean = 101325.0 +prob.meanFlowMag = 0.0 +prob.meanFlowDir = 1 +prob.composition_test = true +prob.Y_CO2_0 = 0.15 # Initial CO2 mass frac contribution +prob.Y_N2_0 = 0.25 # Initial N2 mass frac contribution +prob.Y_AR_0 = 0.6 # Initial AR mass frac contribution +prob.extRhoYCO2 = 145.0 # Source strength for rhoYCO2 +prob.extRhoYCO2_ts = 0.003 # time when source starts + +#---------------------- Temporal CONTROL ------------------------- +peleLM.do_temporals = 1 +peleLM.temporal_int = 1 +peleLM.do_extremas = 1 + +#---------------------------TIME STEPPING------------------------ +amr.stop_time = 0.01 +amr.fixed_dt = 0.001 + +#-------------------------PELE CONTROLS---------------------- +peleLM.v = 1 +peleLM.deltaT_verbose = 1 +peleLM.print_chi_convergence = 1 +peleLM.sdc_iterMax = 2 +peleLM.user_defined_ext_sources = true + +#---------------------------IO CONTROL--------------------------- +amr.check_file = "chk" +amr.check_per = -1 +amr.plot_file = "plt" +amr.plot_per = 0.1 +amr.derive_plot_vars = mass_fractions + +#------------------------- EB SETUP ----------------------------- +eb2.geom_type = "all_regular" \ No newline at end of file diff --git a/Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp b/Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp index 5ae43561..b03e3d66 100644 --- a/Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp +++ b/Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp @@ -26,19 +26,13 @@ prob.T_mean = 300.0 prob.P_mean = 101325.0 prob.meanFlowMag = 4.255 prob.meanFlowDir = 1 +prob.ode_qty_test = true prob.ode_IC = 10.0 prob.ode_xy_lo = -0.01 prob.ode_length = 0.01 prob.ode_height = 0.02 prob.ode_srcstrength = -10.0 -#-----------------INPUTS TO CONSTANT TRANSPORT ------------------ -transport.units = MKS -transport.const_viscosity = 0.0001 -transport.const_bulk_viscosity = 0.0 -transport.const_conductivity = 0.0 -transport.const_diffusivity = 0.0 - #---------------------------TIME STEPPING------------------------ #amr.max_step = 1 amr.stop_time = 0.02 @@ -51,7 +45,6 @@ peleLM.v = 1 peleLM.user_defined_ext_sources = 1 #---------------------------IO CONTROL--------------------------- -#amr.restart = chk_07136 amr.check_file = "chk_" amr.check_per = -1 amr.plot_file = "plt_" diff --git a/Exec/RegTests/EB_ODEQty/eb-odeqty-3d.inp b/Exec/RegTests/EB_ODEQty/eb-odeqty-3d.inp index d839ee73..9c1bb62a 100644 --- a/Exec/RegTests/EB_ODEQty/eb-odeqty-3d.inp +++ b/Exec/RegTests/EB_ODEQty/eb-odeqty-3d.inp @@ -26,19 +26,13 @@ prob.T_mean = 300.0 prob.P_mean = 101325.0 prob.meanFlowMag = 4.255 prob.meanFlowDir = 1 +prob.ode_qty_test = true prob.ode_IC = 10.0 prob.ode_xy_lo = -0.01 prob.ode_length = 0.01 prob.ode_height = 0.02 prob.ode_srcstrength = -10.0 -#-----------------INPUTS TO CONSTANT TRANSPORT ------------------ -transport.units = MKS -transport.const_viscosity = 0.0001 -transport.const_bulk_viscosity = 0.0 -transport.const_conductivity = 0.0 -transport.const_diffusivity = 0.0 - #---------------------------TIME STEPPING------------------------ #amr.max_step = 1 amr.stop_time = 0.02 @@ -51,7 +45,6 @@ peleLM.v = 1 peleLM.user_defined_ext_sources = 1 #---------------------------IO CONTROL--------------------------- -#amr.restart = chk_07136 amr.check_file = "chk_" amr.check_per = -1 amr.plot_file = "plt_" diff --git a/Exec/RegTests/EB_ODEQty/pelelmex_prob.H b/Exec/RegTests/EB_ODEQty/pelelmex_prob.H index 70058e74..49122096 100644 --- a/Exec/RegTests/EB_ODEQty/pelelmex_prob.H +++ b/Exec/RegTests/EB_ODEQty/pelelmex_prob.H @@ -27,8 +27,9 @@ pelelmex_initdata( auto eos = pele::physics::PhysicsType::eos(); amrex::Real massfrac[NUM_SPECIES] = {0.0}; - massfrac[O2_ID] = 0.233; - massfrac[N2_ID] = 0.767; + massfrac[CO2_ID] = prob_parm.Y_CO2_0; + massfrac[N2_ID] = prob_parm.Y_N2_0; + massfrac[AR_ID] = prob_parm.Y_AR_0; switch (prob_parm.meanFlowDir) { case 1: @@ -90,29 +91,31 @@ pelelmex_initdata( #if NUM_ODE > 0 amrex::Real ode_qty[NUM_ODE] = {0.0}; - // Initial state for ode quantities - amrex::Real ode_qty_0[NUM_ODE] = {0.0}; - for (int n = 0; n < NUM_ODE; n++) { - ode_qty_0[n] = (n + 1) * prob_parm.ode_IC; - } - - // Current x,y,z locations - const amrex::Real* prob_lo = geomdata.ProbLo(); - const amrex::Real* dx = geomdata.CellSize(); - const amrex::Real x[AMREX_SPACEDIM] = {AMREX_D_DECL( - prob_lo[0] + static_cast(i + 0.5) * dx[0], - prob_lo[1] + static_cast(j + 0.5) * dx[1], - prob_lo[2] + static_cast(k + 0.5) * dx[2])}; + if (prob_parm.ode_qty_test) { + // Initial state for ode quantities + amrex::Real ode_qty_0[NUM_ODE] = {0.0}; + for (int n = 0; n < NUM_ODE; n++) { + ode_qty_0[n] = (n + 1) * prob_parm.ode_IC; + } - // Create NUM_ODE boxes from (x_strt,y_strt) to (x_end,y_end) - for (int n = 0; n < NUM_ODE; n++) { - amrex::Real x_strt = prob_parm.ode_xy_lo + n * 4. * prob_parm.ode_length; - amrex::Real x_end = x_strt + prob_parm.ode_length; - amrex::Real y_strt = prob_parm.ode_xy_lo; - amrex::Real y_end = y_strt + prob_parm.ode_height; - if ((x[0] > x_strt) && (x[1] > y_strt)) { - if ((x[0] < x_end) && (x[1] < y_end)) { - ode_qty[n] = ode_qty_0[n]; + // Current x,y,z locations + const amrex::Real* prob_lo = geomdata.ProbLo(); + const amrex::Real* dx = geomdata.CellSize(); + const amrex::Real x[AMREX_SPACEDIM] = {AMREX_D_DECL( + prob_lo[0] + static_cast(i + 0.5) * dx[0], + prob_lo[1] + static_cast(j + 0.5) * dx[1], + prob_lo[2] + static_cast(k + 0.5) * dx[2])}; + + // Create NUM_ODE boxes from (x_strt,y_strt) to (x_end,y_end) + for (int n = 0; n < NUM_ODE; n++) { + amrex::Real x_strt = prob_parm.ode_xy_lo + n * 4. * prob_parm.ode_length; + amrex::Real x_end = x_strt + prob_parm.ode_length; + amrex::Real y_strt = prob_parm.ode_xy_lo; + amrex::Real y_end = y_strt + prob_parm.ode_height; + if ((x[0] > x_strt) && (x[1] > y_strt)) { + if ((x[0] < x_end) && (x[1] < y_end)) { + ode_qty[n] = ode_qty_0[n]; + } } } } @@ -139,6 +142,9 @@ bcnormal( pele::physics::PMF::PmfData::DataContainer const* /*pmf_data*/) { amrex::Real massfrac[NUM_SPECIES] = {0.0}; + massfrac[CO2_ID] = prob_parm.Y_CO2_0; + massfrac[N2_ID] = prob_parm.Y_N2_0; + massfrac[AR_ID] = prob_parm.Y_AR_0; auto eos = pele::physics::PhysicsType::eos(); @@ -159,9 +165,6 @@ bcnormal( break; } - massfrac[O2_ID] = 0.333; - massfrac[N2_ID] = 0.667; - s_ext[TEMP] = prob_parm.T_mean; amrex::Real rho_cgs, P_cgs, RhoH_temp; diff --git a/Exec/RegTests/EB_ODEQty/pelelmex_prob.cpp b/Exec/RegTests/EB_ODEQty/pelelmex_prob.cpp index 28f99c2e..478f759c 100644 --- a/Exec/RegTests/EB_ODEQty/pelelmex_prob.cpp +++ b/Exec/RegTests/EB_ODEQty/pelelmex_prob.cpp @@ -10,9 +10,16 @@ PeleLM::readProbParm() // NOLINT(readability-make-member-function-const) pp.query("P_mean", prob_parm->P_mean); pp.query("meanFlowDir", prob_parm->meanFlowDir); pp.query("meanFlowMag", prob_parm->meanFlowMag); + pp.query("ode_qty_test", prob_parm->ode_qty_test); pp.query("ode_IC", prob_parm->ode_IC); pp.query("ode_xy_lo", prob_parm->ode_xy_lo); pp.query("ode_length", prob_parm->ode_length); pp.query("ode_height", prob_parm->ode_height); pp.query("ode_srcstrength", prob_parm->ode_srcstrength); + pp.query("composition_test", prob_parm->composition_test); + pp.query("Y_CO2_0", prob_parm->Y_CO2_0); + pp.query("Y_N2_0", prob_parm->Y_N2_0); + pp.query("Y_AR_0", prob_parm->Y_AR_0); + pp.query("extRhoYCO2", prob_parm->extRhoYCO2); + pp.query("extRhoYCO2_ts", prob_parm->extRhoYCO2_ts); } diff --git a/Exec/RegTests/EB_ODEQty/pelelmex_prob_parm.H b/Exec/RegTests/EB_ODEQty/pelelmex_prob_parm.H index ba51a72e..b9fbbf83 100644 --- a/Exec/RegTests/EB_ODEQty/pelelmex_prob_parm.H +++ b/Exec/RegTests/EB_ODEQty/pelelmex_prob_parm.H @@ -12,11 +12,18 @@ struct ProbParm : amrex::Gpu::Managed amrex::Real T_mean = 298.0_rt; amrex::Real P_mean = 101325.0_rt; amrex::Real meanFlowMag = 0.0; - amrex::Real ode_IC = 10.0; - amrex::Real ode_xy_lo = -0.01; - amrex::Real ode_length = 0.01; - amrex::Real ode_height = 0.020; - amrex::Real ode_srcstrength = -10.0; int meanFlowDir = 1; + bool ode_qty_test = false; + amrex::Real ode_IC = 0.0; + amrex::Real ode_xy_lo = 0.0; + amrex::Real ode_length = 0.0; + amrex::Real ode_height = 0.0; + amrex::Real ode_srcstrength = 0.0; + bool composition_test = false; + amrex::Real Y_CO2_0 = 0.0; + amrex::Real Y_N2_0 = 0.25; + amrex::Real Y_AR_0 = 0.75; + amrex::Real extRhoYCO2 = 0.0; + amrex::Real extRhoYCO2_ts = 0.0; }; #endif diff --git a/Exec/RegTests/EB_ODEQty/plot-comps.py b/Exec/RegTests/EB_ODEQty/plot-comps.py new file mode 100644 index 00000000..5e06ab1f --- /dev/null +++ b/Exec/RegTests/EB_ODEQty/plot-comps.py @@ -0,0 +1,115 @@ +import os +import re +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt + +# File paths +file_dir = os.path.dirname(os.path.abspath(__file__)) +data_file = os.path.join(file_dir, "temporals/tempExtremas") +input_file = os.path.join(file_dir, "composition-test-2d.inp") + +# Load data +col_names = ["time", "max_density", "max_rho.Y(AR)", "max_rho.Y(N2)", "max_rho.Y(CO2)"] +data = pd.read_csv(data_file, usecols=col_names, delimiter=',') +var_names = ["AR", "N2", "CO2"] +time = data['time'] + +# Constants +molar_masses = {"AR": 0.040, "N2": 0.028, "CO2": 0.044} +line_width = 2.5 + +# Parse input file for necessary solution parameters +with open(input_file, 'r') as file: + content = file.read() + soln_params = { + "ts": float(re.search(r"prob\.extRhoYCO2_ts\s*=\s*([\d\.]+)", content).group(1)), + "src_strength": float(re.search(r"prob\.extRhoYCO2\s*=\s*([\d\.]+)", content).group(1)) + } +ts, k = soln_params["ts"], soln_params["src_strength"] + +# Density and species mass data +data_rho = data["max_density"] +data_rhoY = data.iloc[:, 2:].copy() +data_rhoY.columns = var_names + +# Initial values +AR_0, N2_0, CO2_0 = data_rhoY.iloc[0].to_dict().values() +rho_0 = data_rho.iloc[0] + +# Mass fractions and number of moles +data_Y = data_rhoY.div(data_rho, axis=0) +data_Y["Sum(Y_m)"] = data_Y.sum(axis=1) +data_moles = data_rhoY.div([molar_masses[comp] for comp in var_names], axis=1) + +# Exact solutions +exact_rhoY = { + "AR": AR_0, + "N2": N2_0, + "CO2": CO2_0 * (1 - np.heaviside(time - ts, 0)) + + CO2_0 * np.exp(k * (time - ts)) * np.heaviside(time - ts, 0) +} +exact_density = rho_0 + CO2_0 * (np.exp(k * (time - ts)) - 1) * np.heaviside(time - ts, 0) +exact_Y = pd.DataFrame({ + key: exact_rhoY[key] / exact_density for key in var_names}, + index=data.index) +exact_moles = pd.DataFrame({ + key: exact_rhoY[key] / molar_masses[key] for key in var_names} + , index=data.index) + +# Error calculations +max_error_density = np.abs(data_rho - exact_density).max() +max_error_Y = np.abs(data_Y.iloc[:, :-1] - exact_Y).max() +max_error_moles = np.abs(data_moles - exact_moles).max() + +# Print errors +print("\n=====================================================") +print("Max absolute errors:") +print(f"rho: {max_error_density:.16e} kg/m^3") +for species in max_error_Y.index: + print(f"Y({species}): {max_error_Y[species]:.16e}") +for species in max_error_moles.index: + print(f"{species}: {max_error_moles[species]:.16e} moles") +print("=====================================================\n") + +# Plot results +fig, axs = plt.subplots(1, 3, figsize=(16, 4.5)) +xticks = np.linspace(time.min(), time.max(), 5) + +# Plot 1: Mass fractions +for var in var_names: + axs[0].plot(time, data_Y[var], label=f"Y({var})", linewidth=line_width) +axs[0].plot(time, data_Y["Sum(Y_m)"], label="Sum(Y$_m$)", linestyle="--", color="black", linewidth=line_width) +axs[0].set_prop_cycle(None) +for i, var in enumerate(var_names): + axs[0].plot(time, exact_Y[var], 'o', markerfacecolor='none', label='Exact solutions' if i == 0 else None) + +axs[0].set_xticks(xticks) +axs[0].set_xlabel('Time (s)') +axs[0].set_ylabel('Mass Fraction') +axs[0].legend(loc='upper left') +axs[0].grid(True) + +# Plot 2: Density +axs[1].plot(time, data_rho, label="Density", color="black", linewidth=line_width) +axs[1].plot(time, exact_density, 'ko', markerfacecolor='none', label='Exact solution') +axs[1].set_xticks(xticks) +axs[1].set_xlabel('Time (s)') +axs[1].set_ylabel('Density (kg/m$^3$)') +axs[1].legend(loc='upper left') +axs[1].grid(True) + +# Plot 3: Moles +for var in var_names: + axs[2].plot(time, data_moles[var], label=f"Moles of {var}", linewidth=line_width) +axs[2].set_prop_cycle(None) +for i, var in enumerate(var_names): + axs[2].plot(time, exact_moles[var], 'o', markerfacecolor='none', label='Exact solutions' if i == 0 else None) +axs[2].set_xticks(xticks) +axs[2].set_xlabel('Time (s)') +axs[2].set_ylabel('Number of Moles') +axs[2].legend(loc='upper left') +axs[2].grid(True) + +plt.tight_layout() +plt.show() diff --git a/Exec/RegTests/EB_ODEQty/test.py b/Exec/RegTests/EB_ODEQty/test.py new file mode 100644 index 00000000..cbb987cd --- /dev/null +++ b/Exec/RegTests/EB_ODEQty/test.py @@ -0,0 +1,106 @@ +import os +import re +import numpy as np +import pandas as pd +import unittest + +class CompTestCase(unittest.TestCase): + """Test composition of species with external sources""" + + def test_composition(self): + """Are the number of moles, mass fractions, and density correct?""" + + # Molar masses (kg/mol) + molar_masses = {"AR": 0.040, "N2": 0.028, "CO2": 0.044} + + # Load the data + file_dir = os.path.abspath(".") + file_name = os.path.join(file_dir, "temporals/tempExtremas") + col_names = ["time", "max_density", "max_rho.Y(AR)", "max_rho.Y(N2)", "max_rho.Y(CO2)"] + var_names = ["AR", "N2", "CO2"] + data = pd.read_csv(file_name, usecols=col_names, delimiter=',') + print("Raw data:\n", data) + + # Ensure numeric types for all data columns + data = data.apply(pd.to_numeric, errors="coerce") + + # Check for NaN values and raise an error if found + if data.isnull().any().any(): + raise ValueError("Data contains NaN values after conversion.") + + time = data["time"] + + # Parse input file for necessary solution parameters + input_file = os.path.join(file_dir, "composition-test-2d.inp") + with open(input_file, 'r') as file: + content = file.read() + soln_params = { + "ts": float(re.search(r"prob\.extRhoYCO2_ts\s*=\s*([\d\.]+)", content).group(1)), + "src_strength": float(re.search(r"prob\.extRhoYCO2\s*=\s*([\d\.]+)", content).group(1)) + } + ts, k = soln_params["ts"], soln_params["src_strength"] + + # Density and species mass data + data_rho = data["max_density"] + data_rhoY = data.iloc[:, 2:].copy() + data_rhoY.columns = var_names + + # Ensure numeric types for all data columns + data_rho = pd.to_numeric(data_rho, errors="coerce") + data_rhoY = data_rhoY.apply(pd.to_numeric, errors="coerce") + + # Initial values + AR_0, N2_0, CO2_0 = data_rhoY.iloc[0].to_dict().values() + rho_0 = data_rho.iloc[0] + + # Mass fractions and number of moles + data_Y = data_rhoY.div(data_rho, axis=0) + data_Y["Sum(Y_m)"] = data_Y.sum(axis=1) + data_moles = data_rhoY.div([molar_masses[comp] for comp in var_names], axis=1) + + # Exact solutions + exact_rhoY = { + "AR": AR_0, + "N2": N2_0, + "CO2": CO2_0 * (1 - np.heaviside(time - ts, 0)) + + CO2_0 * np.exp(k * (time - ts)) * np.heaviside(time - ts, 0) + } + exact_density = rho_0 + CO2_0 * (np.exp(k * (time - ts)) - 1) * np.heaviside(time - ts, 0) + exact_Y = pd.DataFrame({ + key: exact_rhoY[key] / exact_density for key in var_names}, + index=data.index) + exact_moles = pd.DataFrame({ + key: exact_rhoY[key] / molar_masses[key] for key in var_names}, + index=data.index) + + # Error calculations + max_error_density = np.abs(data_rho - exact_density).max() + max_error_Y = np.abs(data_Y.iloc[:, :-1] - exact_Y).max() + max_error_moles = np.abs(data_moles - exact_moles).max() + + # Expected max errors ["AR", "N2", "CO2"] + expected_errors_moles = np.array([1e-15, 1e-15, 1.0]) + expected_errors_massfrac = np.array([1.032e-02, 4.3e-03, 1.47e-02]) + expected_error_rho = 4e-02 + + # Perform tests using np.testing to make sure errors haven't increased + np.testing.assert_array_less( + np.array([max_error_moles[var] for var in var_names]), + expected_errors_moles, + err_msg="Maximum mole errors exceed specified tolerance." + ) + + np.testing.assert_array_less( + np.array([max_error_Y[var] for var in var_names]), + expected_errors_massfrac, + err_msg="Maximum mass fraction errors exceed specified tolerance." + ) + + np.testing.assert_array_less( + np.array([max_error_density]), + np.array([expected_error_rho]), + err_msg="Maximum density error exceeds specified tolerance." + ) + +if __name__ == "__main__": + unittest.main() diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index d5170821..502fda58 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -275,7 +275,7 @@ PeleLM::getExternalSources( auto* ldata_p_new = getLevelDataPtr(lev, a_timestamp_new); auto& ext_src = m_extSource[lev]; problem_modify_ext_sources( - getTime(lev, a_timestamp_new), m_dt, ldata_p_old->state, + getTime(lev, a_timestamp_old), m_dt, ldata_p_old->state, ldata_p_new->state, ext_src, geom[lev].data(), *prob_parm_d); } } diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index 463e80ef..21fa5e81 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -110,6 +110,14 @@ function(add_test_rev TEST_NAME TEST_EXE_DIR) set_tests_properties(${TEST_NAME} PROPERTIES LABELS "regression;verification;no-ci") endfunction(add_test_rev) +# Regression test with case-specific verification test.py +function(add_test_rt TEST_NAME TEST_EXE_DIR) + setup_test() + set(RUNTIME_OPTIONS "amr.max_step=10 peleLM.do_temporals=1 peleLM.do_extremas=1 peleLM.temporal_int=1 ${RUNTIME_OPTIONS}") + add_test(${TEST_NAME} sh -c "${MPI_COMMANDS} ${CURRENT_TEST_EXE} ${MPIEXEC_POSTFLAGS} ${CURRENT_TEST_BINARY_DIR}/${TEST_NAME}.inp ${RUNTIME_OPTIONS} > ${TEST_NAME}.log ${SAVE_GOLDS_COMMAND} ${FCOMPARE_COMMAND} && nosetests ${CURRENT_TEST_BINARY_DIR}/test.py") + set_tests_properties(${TEST_NAME} PROPERTIES TIMEOUT 18000 PROCESSORS ${PELE_NP} WORKING_DIRECTORY "${CURRENT_TEST_BINARY_DIR}/" LABELS "regression;verification" ATTACHED_FILES_ON_FAIL "${CURRENT_TEST_BINARY_DIR}/${TEST_NAME}.log") +endfunction(add_test_rt) + # Regression tests excluded from CI function(add_test_re TEST_NAME TEST_EXE_DIR) add_test_r(${TEST_NAME} ${TEST_EXE_DIR}) @@ -234,5 +242,8 @@ if(NOT PELE_ENABLE_EB) add_test_r(hit-les-${PELE_DIM}d HITDecay) endif() else() + if(PELE_DIM EQUAL 2) + add_test_rt(composition-test-${PELE_DIM}d EB_ODEQty) + endif() add_test_r(eb-odeqty-${PELE_DIM}d EB_ODEQty) endif()