diff --git a/Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt b/Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt new file mode 100644 index 000000000..38b7706f2 --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt @@ -0,0 +1,8 @@ +set(PELE_PHYSICS_EOS_MODEL Fuego) +set(PELE_PHYSICS_CHEMISTRY_MODEL LiDryer) +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) +set(PELE_PHYSICS_ENABLE_RADIATION OFF) +include(BuildExeAndLib) diff --git a/Exec/RegTests/IsothermalSoretChannel/GNUmakefile b/Exec/RegTests/IsothermalSoretChannel/GNUmakefile new file mode 100644 index 000000000..b509bdcda --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/GNUmakefile @@ -0,0 +1,37 @@ +# AMReX +DIM = 2 +COMP = gnu +PRECISION = DOUBLE +USE_EB = FALSE +USE_HYPRE = FALSE + +# Profiling +PROFILE = FALSE +TINY_PROFILE = FALSE +COMM_PROFILE = FALSE +TRACE_PROFILE = FALSE +MEM_PROFILE = FALSE +USE_GPROF = FALSE + +# Performance +USE_MPI = TRUE +USE_OMP = FALSE +USE_CUDA = FALSE +USE_HIP = FALSE +USE_SYCL = FALSE + +# Debugging +DEBUG = FALSE +FSANITIZER = FALSE +THREAD_SANITIZER = FALSE + +# PeleLMeX +USE_EFIELD = FALSE + +# PelePhysics +Chemistry_Model = LiDryer +Eos_Model = Fuego +Transport_Model = Simple + +PELE_HOME ?= ../../.. +include $(PELE_HOME)/Exec/Make.PeleLMeX diff --git a/Exec/RegTests/IsothermalSoretChannel/README.md b/Exec/RegTests/IsothermalSoretChannel/README.md new file mode 100644 index 000000000..695a820e4 --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/README.md @@ -0,0 +1,3 @@ +## Soret Isothermal wall test + +2D case of periodic channel with isothermal walls and Soret effects enabled. Test to ensure zero species flux at walls. \ No newline at end of file diff --git a/Exec/RegTests/IsothermalSoretChannel/input.inp b/Exec/RegTests/IsothermalSoretChannel/input.inp new file mode 100644 index 000000000..d7360016c --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/input.inp @@ -0,0 +1,91 @@ +#---------------------- DOMAIN DEFINITION ------------------------ +geometry.is_periodic = 1 0 # For each dir, 0: non-perio, 1: periodic +geometry.coord_sys = 0 # 0 => cart, 1 => RZ +geometry.prob_lo = 0.0 0.0 0.0 # x_lo y_lo (z_lo) +geometry.prob_hi = 0.1 0.0125 0.0 # x_hi y_hi (z_hi) + +#---------------------- BC FLAGS --------------------------------- +# Interior, Inflow, Outflow, Symmetry, +# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm +peleLM.lo_bc = Interior NoSlipWallIsotherm # bc in x_lo y_lo (z_lo) +peleLM.hi_bc = Interior NoSlipWallIsotherm # bc in x_hi y_hi (z_hi) + +#---------------------- AMR CONTROL ------------------------------ +amr.n_cell = 256 32 # Level 0 number of cells in each direction +amr.max_level = 1 # maximum level number allowed +amr.ref_ratio = 2 2 2 2 # refinement ratio +amr.regrid_int = 2 # how often to regrid +amr.n_error_buf = 1 1 2 2 # number of buffer cells in error est +amr.grid_eff = 0.7 # what constitutes an efficient grid +amr.blocking_factor = 8 # block factor in grid generation (min box size) +amr.max_grid_size = 32 # max box size + +#---------------------- Problem ---------------------------------- +prob.P_mean = 101325.0 +prob.Vin = 10.0 +prob.Tlow = 400.0 +prob.Thigh = 800.0 +prob.phi = 0.5 +#---------------------- Transport -------------------------------- +transport.use_soret = 1 + +#---------------------- PeleLM CONTROL --------------------------- +peleLM.v = 1 # PeleLMeX verbose +peleLM.use_wbar = 1 # Include Wbar term in species diffusion fluxes +peleLM.sdc_iterMax = 2 # Number of SDC iterations +peleLM.num_init_iter = 3 # Number of initial iterations + +#---------------------- Temporal CONTROL ------------------------- +peleLM.do_temporals = 1 # Turn temporals ON/OFF +peleLM.temporal_int = 10 # Frequency of temporals +peleLM.do_extremas = 1 # Compute state extremas +peleLM.do_mass_balance = 1 # Compute mass balance +peleLM.do_species_balance = 1 # Compute species balance + +#---------------------- Time Stepping CONTROL -------------------- +amr.max_step = 2000 # Maximum number of time steps +amr.stop_time = 4.00 # final simulation physical time [s] +amr.max_wall_time = 1 # Maximum simulation run time [h] +amr.cfl = 0.5 # CFL number for hyperbolic system +amr.dt_shrink = 0.01 # Scale back initial timestep +amr.dt_change_max = 1.1 # Maximum dt increase btw successive steps + +#---------------------- IO CONTROL ------------------------------- +#amr.restart = chk01000 # Restart checkpoint file +amr.check_int = 2000 # Frequency of checkpoint output +amr.plot_int = 20 # Frequency of pltfile output +amr.derive_plot_vars = avg_pressure mag_vort mass_fractions mixture_fraction + +#---------------------- Derived CONTROLS ------------------------- +peleLM.fuel_name = H2 +peleLM.mixtureFraction.format = Cantera +peleLM.mixtureFraction.type = mass +peleLM.mixtureFraction.oxidTank = O2:0.233 N2:0.767 +peleLM.mixtureFraction.fuelTank = H2:1.0 + +#---------------------- Reactor CONTROL -------------------------- +peleLM.chem_integrator = "ReactorNull" +#peleLM.use_typ_vals_chem = 1 # Use species/temp typical values in CVODE +#ode.rtol = 1.0e-6 # Relative tolerance of the chemical solve +#ode.atol = 1.0e-5 # Absolute tolerance factor applied on typical values +#cvode.solve_type = denseAJ_direct # CVODE Linear solve type (for Newton direction) +#cvode.max_order = 4 # CVODE max BDF order. + +#---------------------- Linear solver CONTROL -------------------- +mac_proj.verbose = 0 +nodal_proj.verbose = 0 + +#---------------------- Refinement CONTROL------------------------ +amr.refinement_indicators = highT gradT +amr.highT.max_level = 1 +amr.highT.value_greater = 600 +amr.highT.field_name = temp + +amr.gradT.max_level = 2 +amr.gradT.adjacent_difference_greater = 200 +amr.gradT.field_name = temp + +#---------------------- Debug/HPC CONTROL------------------------- +amrex.fpe_trap_invalid = 1 +amrex.fpe_trap_zero = 1 +amrex.fpe_trap_overflow = 1 diff --git a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H new file mode 100644 index 000000000..33c12a843 --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H @@ -0,0 +1,141 @@ +#ifndef PELELM_PROB_H +#define PELELM_PROB_H + +#include "PeleLMeX_Index.H" +#include "pelelmex_prob_parm.H" +#include "PMFData.H" +#include "PelePhysics.H" + +#include +#include +#include + +AMREX_GPU_HOST_DEVICE +AMREX_FORCE_INLINE +void +set_Y_from_Phi(ProbParm const& prob_parm, amrex::Real Y[]) +{ + amrex::Real phi = prob_parm.phi; + auto eos = pele::physics::PhysicsType::eos(); + amrex::Real Xt[NUM_SPECIES] = {0.0}; + amrex::Real a = 0.0; + switch (prob_parm.fuelID) { +#if defined(H2_ID) + case H2_ID: + a = 0.5; + break; +#endif +#if defined(CH4_ID) + case CH4_ID: + a = 2.0; + break; +#endif + default: + amrex::Abort("fuelID not defined"); + } + + Xt[prob_parm.oxidID] = 1.0 / (1.0 + phi / a + 0.79 / 0.21); + Xt[prob_parm.fuelID] = phi * Xt[prob_parm.oxidID] / a; + Xt[prob_parm.bathID] = 1.0 - Xt[prob_parm.oxidID] - Xt[prob_parm.fuelID]; + + eos.X2Y(Xt, Y); +} + +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +void +pelelmex_initdata( + int i, + int j, + int k, + int /*is_incompressible*/, + amrex::Array4 const& state, + amrex::Array4 const& /*aux*/, + amrex::GeometryData const& geomdata, + ProbParm const& prob_parm, + pele::physics::PMF::PmfData::DataContainer const* /*pmf_data*/) +{ + const amrex::Real* prob_lo = geomdata.ProbLo(); + const amrex::Real* prob_hi = geomdata.ProbHi(); + const amrex::Real* dx = geomdata.CellSize(); + + const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const amrex::Real y_low = prob_lo[1]; + + const amrex::Real Lx = prob_hi[0] - prob_lo[0]; + const amrex::Real Ly = prob_hi[1] - prob_lo[1]; + + auto eos = pele::physics::PhysicsType::eos(); + amrex::Real massfrac[NUM_SPECIES] = {0.0}; + + set_Y_from_Phi(prob_parm, massfrac); + + const amrex::Real T_low = prob_parm.Tlow; + const amrex::Real T_high = prob_parm.Thigh; + + state(i, j, k, TEMP) = (y - y_low) / (Ly) * (T_high - T_low) + T_low; + + state(i, j, k, VELX) = prob_parm.Vin; + state(i, j, k, VELY) = 0.0; + + amrex::Real rho_cgs, P_cgs; + P_cgs = prob_parm.P_mean * 10.0; + + eos.PYT2R(P_cgs, massfrac, state(i, j, k, TEMP), rho_cgs); + state(i, j, k, DENSITY) = rho_cgs * 1.0e3; // CGS -> MKS conversion + + eos.TY2H(state(i, j, k, TEMP), massfrac, state(i, j, k, RHOH)); + state(i, j, k, RHOH) *= + 1.0e-4 * state(i, j, k, DENSITY); // CGS -> MKS conversion + + for (int n = 0; n < NUM_SPECIES; n++) { + state(i, j, k, FIRSTSPEC + n) = massfrac[n] * state(i, j, k, DENSITY); + } +} + +AMREX_GPU_HOST_DEVICE +AMREX_FORCE_INLINE +void +bcnormal( + const amrex::Real x[AMREX_SPACEDIM], + const int /*m_nAux*/, + amrex::Real s_ext[NVAR], + const int /*idir*/, + const int sgn, + const amrex::Real time, + amrex::GeometryData const& /*geomdata*/, + ProbParm const& prob_parm, + pele::physics::PMF::PmfData::DataContainer const* /*pmf_data*/) +{ + auto eos = pele::physics::PhysicsType::eos(); + amrex::Real massfrac[NUM_SPECIES] = {0.0}; + + if (sgn == 1) { + s_ext[TEMP] = prob_parm.Tlow; + } else if (sgn == -1) { + s_ext[TEMP] = prob_parm.Thigh; + } +} + +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +void +zero_visc( + int i, + int j, + int k, + amrex::Array4 const& beta, + amrex::GeometryData const& geomdata, + amrex::Box const& domainBox, + const int dir, + const int beta_comp, + const int nComp) +{ + amrex::ignore_unused( + i, j, k, beta, geomdata, domainBox, dir, beta_comp, nComp); + // We treat species when beta_comp == 0 and nComp == NUM_SPECIES + // otherwise this routine could be called for other face diffusivity (Temp, + // velocity, ...) +} +#endif diff --git a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp new file mode 100644 index 000000000..53c005d5c --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp @@ -0,0 +1,33 @@ +#include +#include + +void +PeleLM::readProbParm() +{ + amrex::ParmParse pp("prob"); + + pp.query("P_mean", PeleLM::prob_parm->P_mean); + pp.query("Vin", PeleLM::prob_parm->Vin); + pp.query("Thigh", PeleLM::prob_parm->Thigh); + pp.query("Tlow", PeleLM::prob_parm->Tlow); + + PeleLM::prob_parm->bathID = N2_ID; + PeleLM::prob_parm->oxidID = O2_ID; + // get the fuel name from the peleLM.fuel_name declaration + amrex::ParmParse pp_pele("peleLM"); + std::string fuelName = ""; + pp_pele.get("fuel_name", fuelName); +#if defined(H2_ID) + if (fuelName == "H2") { + PeleLM::prob_parm->fuelID = H2_ID; + } else +#endif +#if defined(CH4_ID) + if (fuelName == "CH4") { + PeleLM::prob_parm->fuelID = CH4_ID; + } else +#endif + { + amrex::Abort("fuel_name not recognised!"); + } +} diff --git a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H new file mode 100644 index 000000000..6ae9189a8 --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H @@ -0,0 +1,21 @@ +#ifndef PELELM_PROB_PARM_H +#define PELELM_PROB_PARM_H + +#include +#include + +using namespace amrex::literals; + +struct ProbParm +{ + amrex::Real P_mean = 101325.0_rt; + amrex::Real Vin = 10.0; + amrex::Real Thigh = 700.0; + amrex::Real Tlow = 300.0; + amrex::Real phi = 1.0; + + int bathID{-1}; + int fuelID{-1}; + int oxidID{-1}; +}; +#endif diff --git a/Source/Efield/PeleLMeX_EFIonDrift.cpp b/Source/Efield/PeleLMeX_EFIonDrift.cpp index 2e805665c..895a3f6a5 100644 --- a/Source/Efield/PeleLMeX_EFIonDrift.cpp +++ b/Source/Efield/PeleLMeX_EFIonDrift.cpp @@ -41,10 +41,10 @@ PeleLM::ionDriftVelocity(std::unique_ptr& advData) auto bcRecPhiV = fetchBCRecArray(PHIV, 1); getDiffusionOp()->computeGradient( GetVecOfArrOfPtrs(gphiVOld), {}, // don't need the laplacian out - GetVecOfConstPtrs(getPhiVVect(AmrOldTime)), bcRecPhiV[0], do_avgDown); + GetVecOfConstPtrs(getPhiVVect(AmrOldTime)), {}, bcRecPhiV[0], do_avgDown); getDiffusionOp()->computeGradient( GetVecOfArrOfPtrs(gphiVNew), {}, // don't need the laplacian out - GetVecOfConstPtrs(getPhiVVect(AmrNewTime)), bcRecPhiV[0], do_avgDown); + GetVecOfConstPtrs(getPhiVVect(AmrNewTime)), {}, bcRecPhiV[0], do_avgDown); //---------------------------------------------------------------- // TODO : this assumes that all the ions are grouped together at th end ... diff --git a/Source/Efield/PeleLMeX_EFNLSolve.cpp b/Source/Efield/PeleLMeX_EFNLSolve.cpp index 7fffb0187..9c01a8ae9 100644 --- a/Source/Efield/PeleLMeX_EFNLSolve.cpp +++ b/Source/Efield/PeleLMeX_EFNLSolve.cpp @@ -53,7 +53,7 @@ PeleLM::implicitNonLinearSolve( auto bcRecPhiV = fetchBCRecArray(PHIV, 1); getDiffusionOp()->computeGradient( getNLgradPhiVVect(), {}, // don't need the laplacian out - GetVecOfConstPtrs(getPhiVVect(AmrOldTime)), bcRecPhiV[0], do_avgDown); + GetVecOfConstPtrs(getPhiVVect(AmrOldTime)), {}, bcRecPhiV[0], do_avgDown); // Stash away a copy of umac for (int lev = 0; lev <= finest_level; ++lev) { @@ -438,7 +438,7 @@ PeleLM::nonLinearResidual( auto bcRecPhiV = fetchBCRecArray(PHIV, 1); getDiffusionOp()->computeGradient( GetVecOfArrOfPtrs(gradPhiVCur), GetVecOfPtrs(laplacian), - GetVecOfConstPtrs(phiV), bcRecPhiV[0], do_avgDown); + GetVecOfConstPtrs(phiV), {}, bcRecPhiV[0], do_avgDown); // Get nE diffusion term Vector diffnE(finest_level + 1); diff --git a/Source/Efield/PeleLMeX_EFPoisson.cpp b/Source/Efield/PeleLMeX_EFPoisson.cpp index f763f0daa..065a71c4a 100644 --- a/Source/Efield/PeleLMeX_EFPoisson.cpp +++ b/Source/Efield/PeleLMeX_EFPoisson.cpp @@ -47,5 +47,5 @@ PeleLM::poissonSolveEF(const TimeStamp& a_time) // Solve for PhiV getDiffusionOp()->diffuse_scalar( GetVecOfPtrs(getPhiVVect(a_time)), 0, GetVecOfConstPtrs(rhsPoisson), 0, {}, - 0, {}, {}, {}, 0, bcRecPhiV, 1, 1, -eps0 * epsr); + 0, {}, {}, {}, 0, bcRecPhiV, 1, 1, -eps0 * epsr, {}); } diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index c2ccf6fad..8f07e3856 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -553,6 +553,14 @@ public: const amrex::Vector>& a_soretfluxes); + void correctIsothermalBoundary( + const TimeStamp& a_time, + const amrex::Vector& a_spec_boundary, + const amrex::Vector>& + a_wbarfluxes, + const amrex::Vector>& + a_soretfluxes); + /** * \brief Compute the explicit cell-centered diffusion term on all levels * by computing the divergence of the diffusion fluxes at Old or New time, @@ -582,7 +590,8 @@ public: a_spwbarfluxes, amrex::Vector const& a_spec, amrex::Vector const& a_rho, - amrex::Vector const& a_beta); + amrex::Vector const& a_beta, + amrex::Vector const& a_boundary); /** * \brief Add the Soret contribution to the face-centered species diffusion @@ -598,7 +607,6 @@ public: a_spfluxes, const amrex::Vector>& a_spsoretfluxes, - amrex::Vector const& a_spec, amrex::Vector const& a_temp, amrex::Vector const& a_beta); @@ -1893,7 +1901,7 @@ public: int m_use_wbar = 1; int m_use_soret = 0; - + int m_soret_boundary_override = 0; // LES Model bool m_do_les = false; bool m_plot_les = false; diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index ba2ef2a45..5c43ced5e 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -94,8 +94,10 @@ PeleLM::computeDifferentialDiffusionTerms( ? Vector>{} : GetVecOfArrOfPtrs(diffData->wbar_fluxes); Vector> soretFluxVec = - (m_use_soret) != 0 ? GetVecOfArrOfPtrs(diffData->soret_fluxes) - : Vector>{}; + ((is_init != 0) || (m_use_soret == 0)) + ? Vector>{} + : GetVecOfArrOfPtrs(diffData->soret_fluxes); + #ifdef AMREX_USE_EB if (m_isothermalEB != 0) { computeDifferentialDiffusionFluxes( @@ -176,12 +178,12 @@ PeleLM::computeDifferentialDiffusionTerms( #ifdef AMREX_USE_EB fluxDivergenceRD( GetVecOfConstPtrs(getSpeciesVect(a_time)), 0, GetVecOfPtrs(diffData->DT), - 0, GetVecOfArrOfPtrs(diffData->soret_fluxes), 0, {}, 0, NUM_SPECIES, 1, - bcRecSpec_d.dataPtr(), -1.0, m_dt); + 0, GetVecOfArrOfPtrs(diffData->soret_fluxes), 0, {}, 0, NUM_SPECIES, + intensiveFluxes, bcRecSpec_d.dataPtr(), -1.0, m_dt); #else fluxDivergence( GetVecOfPtrs(diffData->DT), 0, GetVecOfArrOfPtrs(diffData->soret_fluxes), - 0, NUM_SPECIES, 1, -1.0); + 0, NUM_SPECIES, intensiveFluxes, -1.0); #endif } @@ -324,6 +326,97 @@ PeleLM::adjustSpeciesFluxes( BL_PROFILE("PeleLM::adjustSpeciesFluxes()"); } +void +PeleLM::correctIsothermalBoundary( + const TimeStamp& a_time, + const Vector& a_spec_boundary, + const Vector>& a_wbarfluxes, + const Vector>& a_soretfluxes) +{ + BL_PROFILE("PeleLMeX::correctIsothermalBoundary()"); + auto bcRecSpec = fetchBCRecArray(FIRSTSPEC, NUM_SPECIES); + bool need_explicit_fluxes = a_soretfluxes.empty(); + + Vector> soretfluxes(finest_level + 1); + if (need_explicit_fluxes) { // need to fill the soret fluxes ourselves + for (int lev = 0; lev <= finest_level; lev++) { + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { + soretfluxes[lev][idim] = new MultiFab( + grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); + soretfluxes[lev][idim]->setVal(0.0); + } + } + addSoretTerm( + soretfluxes, soretfluxes, GetVecOfConstPtrs(getTempVect(a_time)), + GetVecOfConstPtrs(getDiffusivityVect(a_time))); + } else { // have the lagged ones, alias to them + for (int lev = 0; lev <= finest_level; lev++) { + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { + soretfluxes[lev][idim] = new MultiFab( + *a_soretfluxes[lev][idim], amrex::make_alias, 0, NUM_SPECIES); + } + } + } + for (int lev = 0; lev <= finest_level; ++lev) { + auto* ldata_p = getLevelDataPtr(lev, a_time); + // Lets overwrite the boundaries with the fluxes for the inhomogeneous + // Neumann solve + // Get the edge centered diffusivities + MultiFab& ldata_beta_cc = ldata_p->diff_cc; + const Box& domain = geom[lev].Domain(); + int doZeroVisc = 1; + int addTurbContribution = 0; + Array beta_ec = getDiffusivity( + lev, 0, NUM_SPECIES, doZeroVisc, bcRecSpec, ldata_beta_cc, + addTurbContribution); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(ldata_beta_cc, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + const auto bc_lo = m_phys_bc.lo(idim); + const auto bc_hi = m_phys_bc.hi(idim); + const Box& edomain = amrex::surroundingNodes(domain, idim); + const Box& ebx = mfi.nodaltilebox(idim); + auto const& flux_soret = soretfluxes[lev][idim]->const_array(mfi); + auto const& rhoD_ec = beta_ec[idim].const_array(mfi); + auto const& flux_wbar = (m_use_wbar != 0 && !need_explicit_fluxes) + ? a_wbarfluxes[lev][idim]->const_array(mfi) + : rhoD_ec; + auto const& boundary_ar = a_spec_boundary[lev]->array(mfi); + amrex::ParallelFor( + ebx, + [flux_wbar, flux_soret, rhoD_ec, boundary_ar, idim, edomain, bc_lo, + bc_hi, use_wbar = m_use_wbar, + need_explicit_fluxes] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int idx[3] = {i, j, k}; + bool on_lo = (bc_lo == BoundaryCondition::BCNoSlipWallIsotherm || + bc_lo == BoundaryCondition::BCSlipWallIsotherm) && + (idx[idim] <= edomain.smallEnd(idim)); + bool on_hi = (bc_hi == BoundaryCondition::BCNoSlipWallIsotherm || + bc_hi == BoundaryCondition::BCSlipWallIsotherm) && + (idx[idim] >= edomain.bigEnd(idim)); + if (on_lo || on_hi) { + if (on_lo) { // need to move -1 for lo boundary + idx[idim] -= 1; + } + for (int n = 0; n < NUM_SPECIES; n++) { + boundary_ar(idx[0], idx[1], idx[2], n) = flux_soret(i, j, k, n); + // add lagged wbar flux + if (use_wbar != 0 && !need_explicit_fluxes) { + boundary_ar(idx[0], idx[1], idx[2], n) += + flux_wbar(i, j, k, n); + } + boundary_ar(idx[0], idx[1], idx[2], n) /= rhoD_ec(i, j, k, n); + } + } + }); + } + } + } + // TODO: wbar fluxes disabled for this case - boundary system becomes complex +} + void PeleLM::computeDifferentialDiffusionFluxes( const TimeStamp& a_time, @@ -341,6 +434,30 @@ PeleLM::computeDifferentialDiffusionFluxes( BL_PROFILE("PeleLMeX::computeDifferentialDiffusionFluxes()"); //---------------------------------------------------------------- + + Vector spec_boundary; + if (m_soret_boundary_override != 0) { + spec_boundary.resize(finest_level + 1); + // this is the same regardless of lagged or not + for (int lev = 0; lev <= finest_level; ++lev) { + auto* ldata_p = getLevelDataPtr(lev, a_time); + spec_boundary[lev].define( + grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); + + // if we have a mix of Dirichlet and Isothermal walls, we need to give + // Dirichlet boundaries and divide by density since diffuse_scalar doesn't + // touch this boundary MF + + MultiFab::Copy( + spec_boundary[lev], ldata_p->state, FIRSTSPEC, 0, NUM_SPECIES, 1); + for (int n = 0; n < NUM_SPECIES; n++) { + MultiFab::Divide(spec_boundary[lev], ldata_p->state, DENSITY, n, 1, 1); + } + } + // correct the boundary values, pass empty soret & wbar to trigger explicit + // calculation + correctIsothermalBoundary(a_time, GetVecOfPtrs(spec_boundary), {}, {}); + } // Species fluxes // Get the species BCRec auto bcRecSpec = fetchBCRecArray(FIRSTSPEC, NUM_SPECIES); @@ -354,7 +471,7 @@ PeleLM::computeDifferentialDiffusionFluxes( a_fluxes, 0, GetVecOfConstPtrs(getSpeciesVect(a_time)), 0, GetVecOfConstPtrs(getDensityVect(a_time)), GetVecOfConstPtrs(getDiffusivityVect(a_time)), 0, bcRecSpec, - NUM_SPECIES - NUM_IONS, do_avgDown); + NUM_SPECIES - NUM_IONS, do_avgDown, {}); // Ions one by one for (int n = 0; n < NUM_IONS; n++) { auto bcRecIons = fetchBCRecArray(FIRSTSPEC + NUM_SPECIES - NUM_IONS + n, 1); @@ -363,7 +480,7 @@ PeleLM::computeDifferentialDiffusionFluxes( GetVecOfConstPtrs(getSpeciesVect(a_time)), NUM_SPECIES - NUM_IONS + n, GetVecOfConstPtrs(getDensityVect(a_time)), GetVecOfConstPtrs(getDiffusivityVect(a_time)), NUM_SPECIES - NUM_IONS + n, - bcRecIons, 1, do_avgDown); + bcRecIons, 1, do_avgDown, {}); } #else // Get the species diffusion fluxes from the DiffusionOp @@ -374,7 +491,7 @@ PeleLM::computeDifferentialDiffusionFluxes( a_fluxes, 0, GetVecOfConstPtrs(getSpeciesVect(a_time)), 0, GetVecOfConstPtrs(getDensityVect(a_time)), GetVecOfConstPtrs(getDiffusivityVect(a_time)), 0, bcRecSpec, NUM_SPECIES, - do_avgDown); + do_avgDown, GetVecOfConstPtrs(spec_boundary)); #endif // Add the wbar term @@ -384,12 +501,14 @@ PeleLM::computeDifferentialDiffusionFluxes( addWbarTerm( a_fluxes, {}, GetVecOfConstPtrs(getSpeciesVect(a_time)), GetVecOfConstPtrs(getDensityVect(a_time)), - GetVecOfConstPtrs(getDiffusivityVect(a_time))); + GetVecOfConstPtrs(getDiffusivityVect(a_time)), + GetVecOfConstPtrs(spec_boundary)); } else { addWbarTerm( a_fluxes, a_wbarfluxes, GetVecOfConstPtrs(getSpeciesVect(a_time)), GetVecOfConstPtrs(getDensityVect(a_time)), - GetVecOfConstPtrs(getDiffusivityVect(a_time))); + GetVecOfConstPtrs(getDiffusivityVect(a_time)), + GetVecOfConstPtrs(spec_boundary)); } } @@ -398,13 +517,11 @@ PeleLM::computeDifferentialDiffusionFluxes( int need_soret_fluxes = (a_soretfluxes.empty()) ? 0 : 1; if (need_soret_fluxes == 0) { addSoretTerm( - a_fluxes, {}, GetVecOfConstPtrs(getSpeciesVect(a_time)), - GetVecOfConstPtrs(getTempVect(a_time)), + a_fluxes, {}, GetVecOfConstPtrs(getTempVect(a_time)), GetVecOfConstPtrs(getDiffusivityVect(a_time))); } else { addSoretTerm( - a_fluxes, a_soretfluxes, GetVecOfConstPtrs(getSpeciesVect(a_time)), - GetVecOfConstPtrs(getTempVect(a_time)), + a_fluxes, a_soretfluxes, GetVecOfConstPtrs(getTempVect(a_time)), GetVecOfConstPtrs(getDiffusivityVect(a_time))); } } @@ -427,7 +544,6 @@ PeleLM::computeDifferentialDiffusionFluxes( // Set up EB dirichlet value and diffusivity Vector EBvalue(finest_level + 1); Vector EBdiff(finest_level + 1); - ; EBdiff.reserve(finest_level + 1); for (int lev = 0; lev <= finest_level; ++lev) { EBvalue[lev].define( @@ -441,14 +557,14 @@ PeleLM::computeDifferentialDiffusionFluxes( GetVecOfConstPtrs(getTempVect(a_time)), 0, {}, GetVecOfConstPtrs(getDiffusivityVect(a_time)), NUM_SPECIES, GetVecOfConstPtrs(EBvalue), GetVecOfConstPtrs(EBdiff), bcRecTemp, 1, - do_avgDown); + do_avgDown, {}); } else #endif { getDiffusionOp()->computeDiffFluxes( a_fluxes, NUM_SPECIES, GetVecOfConstPtrs(getTempVect(a_time)), 0, {}, GetVecOfConstPtrs(getDiffusivityVect(a_time)), NUM_SPECIES, bcRecTemp, 1, - do_avgDown); + do_avgDown, {}); } // Differential diffusion term: \sum_k ( h_k * \Flux_k ) @@ -477,21 +593,29 @@ PeleLM::addWbarTerm( const Vector>& a_spwbarfluxes, Vector const& a_spec, Vector const& a_rho, - Vector const& a_beta) + Vector const& a_beta, + Vector const& a_boundary) { //------------------------------------------------------------------------ // if a container for wbar fluxes is provided, fill it int need_wbar_fluxes = (a_spwbarfluxes.empty()) ? 0 : 1; - + int have_boundary = (a_boundary.empty()) ? 0 : 1; //------------------------------------------------------------------------ // Compute Wbar on all the levels int nGrow = 1; // Need one ghost cell to compute gradWbar Vector Wbar(finest_level + 1); + Vector Wbar_boundary; + if (have_boundary != 0) { + Wbar_boundary.resize(finest_level + 1); + } auto const* leosparm = eos_parms.device_parm(); for (int lev = 0; lev <= finest_level; ++lev) { - Wbar[lev].define(grids[lev], dmap[lev], 1, nGrow, MFInfo(), Factory(lev)); - + if (have_boundary != 0) { + Wbar_boundary[lev].define( + grids[lev], dmap[lev], 1, nGrow, MFInfo(), Factory(lev)); + } + const Box& domain = geom[lev].Domain(); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -500,10 +624,40 @@ PeleLM::addWbarTerm( auto const& rho_arr = a_rho[lev]->const_array(mfi); auto const& rhoY_arr = a_spec[lev]->const_array(mfi); auto const& Wbar_arr = Wbar[lev].array(mfi); + auto const& gradY_arr = + (have_boundary != 0) ? a_boundary[lev]->const_array(mfi) : Wbar_arr; + auto const& Wbar_boundary_arr = + (have_boundary != 0) ? Wbar_boundary[lev].array(mfi) : Wbar_arr; + amrex::ParallelFor( - gbx, [rho_arr, rhoY_arr, Wbar_arr, + + gbx, [rho_arr, rhoY_arr, Wbar_arr, gradY_arr, Wbar_boundary_arr, domain, + have_boundary, phys_bc = m_phys_bc, leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { getMwmixGivenRY(i, j, k, rho_arr, rhoY_arr, Wbar_arr, leosparm); + if (have_boundary != 0) { // need to impose gradWbar on boundary for + // computeGradient + // for dirichlet boundaries, we'll overwrite inhomog neumann ones + // NOTE: for now, this is skipped since wbar disabled for + // isothermal/soret + Wbar_boundary_arr(i, j, k) = Wbar_arr(i, j, k); + int idx[3] = {i, j, k}; + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { + const auto bc_lo = phys_bc.lo(idim); + const auto bc_hi = phys_bc.hi(idim); + bool on_lo = (bc_lo == BoundaryCondition::BCNoSlipWallIsotherm || + bc_lo == BoundaryCondition::BCSlipWallIsotherm) && + (idx[idim] < domain.smallEnd(idim)); + bool on_hi = (bc_hi == BoundaryCondition::BCNoSlipWallIsotherm || + bc_hi == BoundaryCondition::BCSlipWallIsotherm) && + (idx[idim] > domain.bigEnd(idim)); + + if (on_lo || on_hi) { + getGradMwmixGivengradYMwmix( + i, j, k, gradY_arr, Wbar_arr, Wbar_boundary_arr, leosparm); + } + } + } }); } } @@ -528,7 +682,8 @@ PeleLM::addWbarTerm( } getDiffusionOp()->computeGradient( GetVecOfArrOfPtrs(gradWbar), {}, // Don't need the laplacian out - GetVecOfConstPtrs(Wbar), bcRecSpec[0], do_avgDown); + GetVecOfConstPtrs(Wbar), GetVecOfConstPtrs(Wbar_boundary), bcRecSpec[0], + do_avgDown); //------------------------------------------------------------------------ // add Wbar term to species fluxes @@ -583,8 +738,9 @@ PeleLM::addWbarTerm( ? a_spwbarfluxes[lev][idim]->array(mfi) : a_spfluxes[lev][idim]->array(mfi); // Dummy unused Array4 - // Wbar flux is : - \rho Y_m / \overline{W} * D_m * \nabla - // \overline{W} with beta_m = \rho * D_m below + // Wbar flux is : - \rho Y_m / W_k * D_m * \nabla + // \overline{W} with beta_m = \rho * D_m * overline(W) / W_k below + // need to divide by \overline(W) amrex::ParallelFor( ebx, [need_wbar_fluxes, gradWbar_ar, beta_ar, rhoY, spFlux_ar, spwbarFlux_ar, @@ -625,7 +781,6 @@ void PeleLM::addSoretTerm( const Vector>& a_spfluxes, const Vector>& a_spsoretfluxes, - Vector const& a_spec, Vector const& a_temp, Vector const& a_beta) { @@ -654,7 +809,7 @@ PeleLM::addSoretTerm( } getDiffusionOp()->computeGradient( GetVecOfArrOfPtrs(gradT), {}, // Don't need the laplacian out - a_temp, bcRecTemp[0], do_avgDown); + a_temp, {}, bcRecTemp[0], do_avgDown); //------------------------------------------------------------------------ // add Soret term to species fluxes @@ -672,40 +827,17 @@ PeleLM::addSoretTerm( #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif { - FArrayBox rhoY_ed; FArrayBox T_ed; for (MFIter mfi(*a_beta[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { // Get edge centered rhoYs const Box ebx = mfi.nodaltilebox(idim); - rhoY_ed.resize(ebx, NUM_SPECIES); - Elixir rhoY_el = rhoY_ed.elixir(); - const Box& edomain = amrex::surroundingNodes(domain, idim); - auto const& rhoY_arr = a_spec[lev]->const_array(mfi); - const auto& rhoYed_arr = rhoY_ed.array(0); - const auto bc_lo_spec = bcRecSpec[0].lo(idim); - const auto bc_hi_spec = bcRecSpec[0].hi(idim); - amrex::ParallelFor( - ebx, [idim, bc_lo_spec, bc_hi_spec, use_harmonic_avg, rhoY_arr, - rhoYed_arr, - edomain] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - int idx[3] = {i, j, k}; - bool on_lo = - ((bc_lo_spec == amrex::BCType::ext_dir) && - (idx[idim] <= edomain.smallEnd(idim))); - bool on_hi = - ((bc_hi_spec == amrex::BCType::ext_dir) && - (idx[idim] >= edomain.bigEnd(idim))); - cen2edg_cpp( - i, j, k, idim, NUM_SPECIES, use_harmonic_avg, on_lo, on_hi, - rhoY_arr, rhoYed_arr); - }); // Get edge centered temps T_ed.resize(ebx, 1); - Elixir T_el = T_ed.elixir(); // point of this? + Elixir T_el = T_ed.elixir(); auto const& T_arr = a_temp[lev]->const_array(mfi); const auto& Ted_arr = T_ed.array(0); @@ -736,8 +868,8 @@ PeleLM::addSoretTerm( ? a_spsoretfluxes[lev][idim]->array(mfi) : a_spfluxes[lev][idim]->array(mfi); // Dummy unused Array4 - // Soret flux is : - \rho * Y theta_m * \nabla T / T - // with beta_m = \rho (* Y?) * theta_m below + // Soret flux is : - rho * D_m * chi_m * \nabla T / T + // with beta_m = rho * D_m * chi_m below amrex::ParallelFor( ebx, [need_soret_fluxes, gradT_ar, beta_ar, T, spFlux_ar, @@ -919,6 +1051,32 @@ PeleLM::differentialDiffusionUpdate( // Get the species BCRec auto bcRecSpec = fetchBCRecArray(FIRSTSPEC, NUM_SPECIES); + Vector spec_boundary; + if (m_soret_boundary_override != 0) { + spec_boundary.resize(finest_level + 1); + // this is the same regardless of lagged or not + for (int lev = 0; lev <= finest_level; ++lev) { + auto* ldata_p = getLevelDataPtr(lev, AmrNewTime); + spec_boundary[lev].define( + grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); + + // if we have a mix of Dirichlet and Isothermal walls, we need to give + // Dirichlet boundaries and divide by density since diffuse_scalar doesn't + // touch this boundary MF + + MultiFab::Copy( + spec_boundary[lev], ldata_p->state, FIRSTSPEC, 0, NUM_SPECIES, 1); + for (int n = 0; n < NUM_SPECIES; n++) { + MultiFab::Divide(spec_boundary[lev], ldata_p->state, DENSITY, n, 1, 1); + } + } + // correct boundary, we have lagged fluxes so lets use them + correctIsothermalBoundary( + AmrNewTime, GetVecOfPtrs(spec_boundary), + GetVecOfArrOfPtrs(diffData->wbar_fluxes), + GetVecOfArrOfPtrs(diffData->soret_fluxes)); + } + #ifdef PELE_USE_EFIELD // Solve for \widetilda{rhoY^{np1,kp1}} // -> return the uncorrected fluxes^{np1,kp1} @@ -933,7 +1091,7 @@ PeleLM::differentialDiffusionUpdate( GetVecOfConstPtrs( getDensityVect(AmrNewTime)), // this triggers proper scaling by density GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), 0, bcRecSpec, - NUM_SPECIES - NUM_IONS, 0, m_dt); + NUM_SPECIES - NUM_IONS, 0, m_dt, {}); // Ions one by one for (int n = 0; n < NUM_IONS; n++) { auto bcRecIons = fetchBCRecArray(FIRSTSPEC + NUM_SPECIES - NUM_IONS + n, 1); @@ -946,7 +1104,7 @@ PeleLM::differentialDiffusionUpdate( GetVecOfConstPtrs( getDensityVect(AmrNewTime)), // this triggers proper scaling by density GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), 0, bcRecIons, 1, 0, - m_dt); + m_dt, {}); } #else // Solve for \widetilda{rhoY^{np1,kp1}} @@ -962,7 +1120,7 @@ PeleLM::differentialDiffusionUpdate( GetVecOfConstPtrs( getDensityVect(AmrNewTime)), // this triggers proper scaling by density GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), 0, bcRecSpec, - NUM_SPECIES, 0, m_dt); + NUM_SPECIES, 0, m_dt, GetVecOfConstPtrs(spec_boundary)); #endif // Add lagged Wbar term @@ -1104,7 +1262,7 @@ PeleLM::differentialDiffusionUpdate( GetVecOfConstPtrs(getTempVect(AmrNewTime)), 0, {}, GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), NUM_SPECIES, GetVecOfConstPtrs(EBvalue), GetVecOfConstPtrs(EBdiff), bcRecTemp, 1, - do_avgDown); + do_avgDown, {}); } else #endif { @@ -1112,7 +1270,7 @@ PeleLM::differentialDiffusionUpdate( GetVecOfArrOfPtrs(fluxes), NUM_SPECIES, GetVecOfConstPtrs(getTempVect(AmrNewTime)), 0, {}, GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), NUM_SPECIES, bcRecTemp, - 1, do_avgDown); + 1, do_avgDown, {}); } // Differential diffusion term: \sum_k ( h_k * \Flux_k ) @@ -1193,7 +1351,7 @@ PeleLM::differentialDiffusionUpdate( GetVecOfConstPtrs(rhs), 0, GetVecOfArrOfPtrs(fluxes), NUM_SPECIES, GetVecOfConstPtrs(RhoCp), {}, GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), NUM_SPECIES, - GetVecOfConstPtrs(EBdiff), 0, bcRecTemp, 1, 0, m_dt); + GetVecOfConstPtrs(EBdiff), 0, bcRecTemp, 1, 0, m_dt, {}); } else #endif { @@ -1201,7 +1359,7 @@ PeleLM::differentialDiffusionUpdate( GetVecOfPtrs(getTempVect(AmrNewTime)), 0, GetVecOfConstPtrs(rhs), 0, GetVecOfArrOfPtrs(fluxes), NUM_SPECIES, GetVecOfConstPtrs(RhoCp), {}, GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), NUM_SPECIES, - bcRecTemp, 1, 0, m_dt); + bcRecTemp, 1, 0, m_dt, {}); } // Post deltaT iteration linear solve @@ -1349,14 +1507,14 @@ PeleLM::deltaTIter_update( GetVecOfConstPtrs(getTempVect(AmrNewTime)), 0, {}, GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), NUM_SPECIES, GetVecOfConstPtrs(EBvalue), GetVecOfConstPtrs(EBdiff), bcRecTemp, 1, - do_avgDown); + do_avgDown, {}); } else #endif { getDiffusionOp()->computeDiffFluxes( a_fluxes, NUM_SPECIES, GetVecOfConstPtrs(getTempVect(AmrNewTime)), 0, {}, GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), NUM_SPECIES, bcRecTemp, - 1, do_avgDown); + 1, do_avgDown, {}); } // Differential diffusion term: \sum_k ( h_k * \Flux_k ) @@ -1535,6 +1693,19 @@ PeleLM::getDiffusionLinOpBC(Orientation::Side a_side, const BCRec& a_bc) amrexbc == amrex::BCType::hoextrap || amrexbc == amrex::BCType::reflect_even) { r[idim] = LinOpBCType::Neumann; + if ( + m_use_soret != 0 && + amrexbc == amrex::BCType::foextrap) { // should just catch + // density/species/forcing for + // isothermal walls + auto phys_bc = (a_side == Orientation::low) ? m_phys_bc.lo(idim) + : m_phys_bc.hi(idim); + if ( + phys_bc == BoundaryCondition::BCNoSlipWallIsotherm || + phys_bc == BoundaryCondition::BCSlipWallIsotherm) { + r[idim] = LinOpBCType::inhomogNeumann; + } + } } else if (amrexbc == amrex::BCType::reflect_odd) { r[idim] = LinOpBCType::reflect_odd; } else { diff --git a/Source/PeleLMeX_DiffusionOp.H b/Source/PeleLMeX_DiffusionOp.H index abeb9f494..d4705b37d 100644 --- a/Source/PeleLMeX_DiffusionOp.H +++ b/Source/PeleLMeX_DiffusionOp.H @@ -36,7 +36,8 @@ public: amrex::Vector a_bcrec, int ncomp, int isPoissonSolve, - amrex::Real dt); + amrex::Real dt, + amrex::Vector const& a_boundary); #ifdef AMREX_USE_EB void diffuse_scalar( @@ -58,7 +59,8 @@ public: amrex::Vector a_bcrec, int ncomp, int isPoissonSolve, - amrex::Real dt); + amrex::Real dt, + amrex::Vector const& a_boundary); #endif void computeDiffLap( @@ -82,7 +84,8 @@ public: int bcoeff_comp, amrex::Vector a_bcrec, int ncomp, - int do_avgDown); + int do_avgDown, + amrex::Vector const& a_boundary); #ifdef AMREX_USE_EB void computeDiffFluxes( @@ -100,13 +103,15 @@ public: amrex::Vector const& a_EBbcoeff, amrex::Vector a_bcrec, int ncomp, - int do_avgDown); + int do_avgDown, + amrex::Vector const& a_boundary); #endif void computeGradient( const amrex::Vector>& a_grad, const amrex::Vector& a_laps, const amrex::Vector& a_phi, + const amrex::Vector& a_boundary, const amrex::BCRec& a_bcrec, int do_avgDown) const; diff --git a/Source/PeleLMeX_DiffusionOp.cpp b/Source/PeleLMeX_DiffusionOp.cpp index e931923b9..9ac11506d 100644 --- a/Source/PeleLMeX_DiffusionOp.cpp +++ b/Source/PeleLMeX_DiffusionOp.cpp @@ -106,7 +106,8 @@ DiffusionOp::diffuse_scalar( Vector a_bcrec, int ncomp, int isPoissonSolve, - Real a_dt) + Real a_dt, + Vector const& a_boundary) { BL_PROFILE("DiffusionOp::diffuse_scalar()"); @@ -116,6 +117,7 @@ DiffusionOp::diffuse_scalar( int have_fluxes = (a_flux.empty()) ? 0 : 1; int have_acoeff = (a_acoeff.empty()) ? 0 : 1; int have_bcoeff = (a_bcoeff.empty()) ? 0 : 1; + int have_boundary = (a_boundary.empty()) ? 0 : 1; //---------------------------------------------------------------- // Checks @@ -193,6 +195,7 @@ DiffusionOp::diffuse_scalar( Vector> fluxes(finest_level + 1); Vector component; Vector rhs; + Vector boundary; // Allow for component specific LinOp BC m_scal_solve_op->setDomainBC( @@ -229,7 +232,13 @@ DiffusionOp::diffuse_scalar( component.emplace_back(phi[lev], amrex::make_alias, comp, m_ncomp); rhs.emplace_back( *a_rhs[lev], amrex::make_alias, rhs_comp + comp, m_ncomp); - m_scal_solve_op->setLevelBC(lev, &component[lev]); + if (have_boundary != 0) { + boundary.emplace_back( + *a_boundary[lev], amrex::make_alias, comp, m_ncomp); + } else { + boundary.emplace_back(phi[lev], amrex::make_alias, comp, m_ncomp); + } + m_scal_solve_op->setLevelBC(lev, &boundary[lev]); } // Setup linear solver @@ -317,7 +326,8 @@ DiffusionOp::diffuse_scalar( Vector a_bcrec, int ncomp, int isPoissonSolve, - Real a_dt) + Real a_dt, + Vector const& a_boundary) { BL_PROFILE("DiffusionOp::diffuse_scalar()"); @@ -327,6 +337,7 @@ DiffusionOp::diffuse_scalar( int have_fluxes = (a_flux.empty()) ? 0 : 1; int have_acoeff = (a_acoeff.empty()) ? 0 : 1; int have_bcoeff = (a_bcoeff.empty()) ? 0 : 1; + int have_boundary = (a_boundary.empty()) ? 0 : 1; //---------------------------------------------------------------- // Checks @@ -404,6 +415,7 @@ DiffusionOp::diffuse_scalar( Vector> fluxes(finest_level + 1); Vector component; Vector rhs; + Vector boundary; // Allow for component specific LinOp BC m_scal_solve_op->setDomainBC( @@ -435,7 +447,13 @@ DiffusionOp::diffuse_scalar( component.emplace_back(phi[lev], amrex::make_alias, comp, m_ncomp); rhs.emplace_back( *a_rhs[lev], amrex::make_alias, rhs_comp + comp, m_ncomp); - m_scal_solve_op->setLevelBC(lev, &component[lev]); + if (have_boundary != 0) { + boundary.emplace_back( + *a_boundary[lev], amrex::make_alias, comp, m_ncomp); + } else { + boundary.emplace_back(phi[lev], amrex::make_alias, comp, m_ncomp); + } + m_scal_solve_op->setLevelBC(lev, &boundary[lev]); m_scal_solve_op->setEBDirichlet(lev, *a_phiEB[lev], *a_bcoeffEB[lev]); } @@ -585,7 +603,8 @@ DiffusionOp::computeDiffFluxes( int bcoeff_comp, Vector a_bcrec, int ncomp, - int do_avgDown) + int do_avgDown, + Vector const& a_boundary) { BL_PROFILE("DiffusionOp::computeDiffFluxes()"); @@ -600,6 +619,7 @@ DiffusionOp::computeDiffFluxes( int finest_level = m_pelelm->finestLevel(); int have_density = (a_density.empty()) ? 0 : 1; + int have_boundary = (a_boundary.empty()) ? 0 : 1; // Duplicate phi since it is modified by the LinOp // and if have_density -> divide by density @@ -648,6 +668,7 @@ DiffusionOp::computeDiffFluxes( Vector> fluxes(finest_level + 1); Vector component; Vector laps; + Vector boundary; // Allow for component specific LinOp BC m_scal_apply_op->setDomainBC( @@ -660,6 +681,13 @@ DiffusionOp::computeDiffFluxes( *a_flux[lev][idim], amrex::make_alias, flux_comp + comp, m_ncomp); } component.emplace_back(phi[lev], amrex::make_alias, comp, m_ncomp); + if (have_boundary != 0) { + boundary.emplace_back( + *a_boundary[lev], amrex::make_alias, comp, m_ncomp); + } else { + boundary.emplace_back(phi[lev], amrex::make_alias, comp, m_ncomp); + } + int doZeroVisc = 1; int addTurbContrib = 1; Vector subBCRec = { @@ -676,7 +704,7 @@ DiffusionOp::computeDiffFluxes( #else m_scal_apply_op->setBCoeffs(lev, GetArrOfConstPtrs(bcoeff_ec)); #endif - m_scal_apply_op->setLevelBC(lev, &component[lev]); + m_scal_apply_op->setLevelBC(lev, &boundary[lev]); } MLMG mlmg(*m_scal_apply_op); @@ -716,7 +744,8 @@ DiffusionOp::computeDiffFluxes( Vector const& a_EBbcoeff, Vector a_bcrec, int ncomp, - int do_avgDown) + int do_avgDown, + Vector const& a_boundary) { BL_PROFILE("DiffusionOp::computeDiffFluxes()"); @@ -731,6 +760,7 @@ DiffusionOp::computeDiffFluxes( int finest_level = m_pelelm->finestLevel(); int have_density = (a_density.empty()) ? 0 : 1; + int have_boundary = (a_boundary.empty()) ? 0 : 1; // Duplicate phi since it is modified by the LinOp // and if have_density -> divide by density @@ -781,6 +811,7 @@ DiffusionOp::computeDiffFluxes( Vector> ebfluxes; Vector component; Vector laps; + Vector boundary; // Allow for component specific LinOp BC m_scal_apply_op->setDomainBC( @@ -795,6 +826,12 @@ DiffusionOp::computeDiffFluxes( ebfluxes.push_back(std::make_unique( *a_EBflux[lev], amrex::make_alias, ebflux_comp + comp, m_ncomp)); component.emplace_back(phi[lev], amrex::make_alias, comp, m_ncomp); + if (have_boundary != 0) { + boundary.emplace_back( + *a_boundary[lev], amrex::make_alias, comp, m_ncomp); + } else { + boundary.emplace_back(phi[lev], amrex::make_alias, comp, m_ncomp); + } int doZeroVisc = 1; Vector subBCRec = { a_bcrec.begin() + comp, a_bcrec.begin() + comp + m_ncomp}; @@ -805,7 +842,7 @@ DiffusionOp::computeDiffFluxes( MFInfo(), a_phi[lev]->Factory()); m_scal_apply_op->setBCoeffs( lev, GetArrOfConstPtrs(bcoeff_ec), MLMG::Location::FaceCentroid); - m_scal_apply_op->setLevelBC(lev, &component[lev]); + m_scal_apply_op->setLevelBC(lev, &boundary[lev]); m_scal_apply_op->setEBDirichlet(lev, *a_EBvalue[lev], *a_EBbcoeff[lev]); } @@ -829,6 +866,7 @@ DiffusionOp::computeGradient( const Vector>& a_grad, const Vector& a_laps, const Vector& a_phi, + const Vector& a_boundary, const BCRec& a_bcrec, int do_avgDown) const { @@ -836,7 +874,6 @@ DiffusionOp::computeGradient( // Do I need the Laplacian out ? int need_laplacian = (a_laps.empty()) ? 0 : 1; - // Force updating the operator for (int lev = 0; lev <= m_pelelm->finestLevel(); ++lev) { m_gradient_op->setBCoeffs(lev, -1.0); @@ -847,6 +884,7 @@ DiffusionOp::computeGradient( AMREX_ASSERT(a_phi[0]->nGrow() >= 1); int finest_level = m_pelelm->finestLevel(); + int have_boundary = (a_boundary.empty()) ? 0 : 1; // Set domainBCs m_gradient_op->setDomainBC( @@ -856,13 +894,25 @@ DiffusionOp::computeGradient( // Duplicate phi since it is modified by the LinOp // and setup level BCs Vector phi(finest_level + 1); + Vector boundary(finest_level + 1); Vector laps; for (int lev = 0; lev <= finest_level; ++lev) { phi[lev].define( a_phi[lev]->boxArray(), a_phi[lev]->DistributionMap(), 1, 1, MFInfo(), a_phi[lev]->Factory()); + boundary[lev].define( + a_phi[lev]->boxArray(), a_phi[lev]->DistributionMap(), 1, 1, MFInfo(), + a_phi[lev]->Factory()); + MultiFab::Copy(phi[lev], *a_phi[lev], 0, 0, 1, 1); - m_gradient_op->setLevelBC(lev, &phi[lev]); + + if (have_boundary != 0) { + MultiFab::Copy(boundary[lev], *a_boundary[lev], 0, 0, 1, 1); + } else { + MultiFab::Copy(boundary[lev], *a_phi[lev], 0, 0, 1, 1); + } + + m_gradient_op->setLevelBC(lev, &boundary[lev]); if (need_laplacian != 0) { laps.emplace_back(*a_laps[lev], amrex::make_alias, 0, 1); } else { diff --git a/Source/PeleLMeX_K.H b/Source/PeleLMeX_K.H index 13f4a1914..3ffbfa0ac 100644 --- a/Source/PeleLMeX_K.H +++ b/Source/PeleLMeX_K.H @@ -562,6 +562,33 @@ getMwmixGivenRY( Mwmix(i, j, k) = Mwmix(i, j, k) * 0.001_rt; // CGS -> MKS conversion } +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +void +getGradMwmixGivengradYMwmix( + int i, + int j, + int k, + amrex::Array4 const& gradY, + amrex::Array4 const& Mwmix, + amrex::Array4 const& gradMwmix, + pele::physics::eos::EosParm const* + eosparm) noexcept +{ + using namespace amrex::literals; + + auto eos = pele::physics::PhysicsType::eos(eosparm); + amrex::Real imw[NUM_SPECIES] = {0.0_rt}; + eos.inv_molecular_weight(imw); + gradMwmix(i, j, k) = 0.0; + for (int n = 0; n < NUM_SPECIES; n++) { + imw[n] *= 1000.0; // CGS -> MKS + gradMwmix(i, j, k) += gradY(i, j, k, n) * imw[n]; + } + // gradWbar = -Wbar^2 * sum(grad Y_k / W_k) + gradMwmix(i, j, k) *= -Mwmix(i, j, k) * Mwmix(i, j, k); +} + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 9240a0cb0..292030ced 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -105,6 +105,12 @@ PeleLM::Setup() amrex::Print() << " Using mixture-averaged transport with Soret effects" << std::endl; + if (m_soret_boundary_override != 0) { + amrex::Print() + << " Imposing inhomogeneous Neumann conditions " + "for species on isothermal walls. WARNING: use_wbar disabled." + << std::endl; + } } } else { if (m_fixed_Le != 0) { @@ -417,6 +423,23 @@ PeleLM::readParameters() ParmParse pptrans("transport"); pptrans.query("use_soret", m_use_soret); pp.query("use_wbar", m_use_wbar); + if (m_use_soret != 0) { + bool isothermal = false; + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { + isothermal |= + (m_phys_bc.lo(idim) == BoundaryCondition::BCSlipWallIsotherm || + m_phys_bc.lo(idim) == BoundaryCondition::BCNoSlipWallIsotherm || + m_phys_bc.hi(idim) == BoundaryCondition::BCSlipWallIsotherm || + m_phys_bc.hi(idim) == BoundaryCondition::BCNoSlipWallIsotherm); + } + if (isothermal) { + m_soret_boundary_override = 1; + m_use_wbar = 0; +#if PELE_USE_EFIELD + amrex::Abort("Isothermal walls with Soret incompatible with Efield"); +#endif + } + } pp.query("unity_Le", m_unity_Le); pp.query("fixed_Le", m_fixed_Le); pp.query("fixed_Pr", m_fixed_Pr); diff --git a/Submodules/AMReX-Hydro b/Submodules/AMReX-Hydro index 3788361fe..13141b537 160000 --- a/Submodules/AMReX-Hydro +++ b/Submodules/AMReX-Hydro @@ -1 +1 @@ -Subproject commit 3788361fec234e29e6d4e8ba95010e68d9c426e7 +Subproject commit 13141b5373d83ef0f30b057976c0c7ac82291f9d