From 3c847ae19864fa9c641cc267d42d86e12d8912bd Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 21 Jun 2024 15:08:37 +0200 Subject: [PATCH 01/34] Soret Isothermal wall fix --- Source/PeleLMeX.H | 2 +- Source/PeleLMeX_Diffusion.cpp | 95 +++++++++++++++++++++++++++++++-- Source/PeleLMeX_DiffusionOp.H | 3 +- Source/PeleLMeX_DiffusionOp.cpp | 12 ++++- Source/PeleLMeX_Init.cpp | 3 +- Source/PeleLMeX_Regrid.cpp | 6 +-- Source/PeleLMeX_Setup.cpp | 18 +++++++ 7 files changed, 126 insertions(+), 13 deletions(-) diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 47fcb7750..db69ffc4a 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -1856,7 +1856,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 18c341f58..2927b10ba 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -176,12 +176,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, + 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 } @@ -903,6 +903,82 @@ PeleLM::differentialDiffusionUpdate( // Get the species BCRec auto bcRecSpec = fetchBCRecArray(FIRSTSPEC, NUM_SPECIES); + Vector spec_boundary(finest_level + 1); + for (int lev = 0; lev <= finest_level; ++lev) { + auto* ldata_p = getLevelDataPtr(lev, AmrNewTime); + spec_boundary[lev] = new MultiFab( + ldata_p->state.boxArray(), ldata_p->state.DistributionMap(), NUM_SPECIES, + 1, MFInfo(), ldata_p->state.Factory()); + // copy the state into the boundaries for Dirichlet boundaries + MultiFab::Copy( + *spec_boundary[lev], ldata_p->state, FIRSTSPEC, 0, NUM_SPECIES, 1); + // to remain consistent with diffusion solve, we need to divide by density + // here if we have mix of Dirichlet and Isothermal conditions, we need to do + // work here, rather than letting nullptr do the work + for (int n = 0; n < NUM_SPECIES; n++) { + MultiFab::Divide(*spec_boundary[lev], ldata_p->state, DENSITY, n, 1, 1); + } + // no Isothermal conditions, just carry Dirichlet conditions through to + // diffusion solve + if (m_soret_boundary_override == 0) { + break; + } + // otherwise, lets overwrite the boundaries with the fluxes + 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(*spec_boundary[lev], 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 = + diffData->soret_fluxes[lev][idim].const_array(mfi); + auto const& flux_wbar = + diffData->wbar_fluxes[lev][idim].const_array(mfi); + auto const& rhoD_ec = beta_ec[idim].const_array(mfi); + auto const& boundary_ar = spec_boundary[lev]->array(mfi); + const int use_wbar = m_use_wbar; + amrex::ParallelFor( + ebx, + [flux_wbar, flux_soret, rhoD_ec, boundary_ar, idim, edomain, bc_lo, + bc_hi, use_wbar] 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) { + 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); + } + } + }); + } + } + } + #ifdef PELE_USE_EFIELD // Solve for \widetilda{rhoY^{np1,kp1}} // -> return the uncorrected fluxes^{np1,kp1} @@ -946,7 +1022,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, spec_boundary); #endif // Add lagged Wbar term @@ -1517,6 +1593,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..a365b01ab 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 a_boundary = {nullptr}); #ifdef AMREX_USE_EB void diffuse_scalar( diff --git a/Source/PeleLMeX_DiffusionOp.cpp b/Source/PeleLMeX_DiffusionOp.cpp index e931923b9..2309bbd92 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 a_boundary) { BL_PROFILE("DiffusionOp::diffuse_scalar()"); @@ -193,6 +194,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 +231,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 (a_boundary[lev] != nullptr) { + 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 diff --git a/Source/PeleLMeX_Init.cpp b/Source/PeleLMeX_Init.cpp index 239b6c30c..fbe0b14de 100644 --- a/Source/PeleLMeX_Init.cpp +++ b/Source/PeleLMeX_Init.cpp @@ -38,8 +38,7 @@ PeleLM::MakeNewLevelFromScratch( auto const dx = geom[lev].CellSizeArray(); Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); amrex::Print() << " with " << ba.numPts() << " cells," << ba.size() - << " boxes," - << " over " + << " boxes," << " over " << static_cast(ba.numPts()) * vol / geom[0].ProbSize() * 100 << "% of the domain \n"; diff --git a/Source/PeleLMeX_Regrid.cpp b/Source/PeleLMeX_Regrid.cpp index c4d3fb2ba..b61a5ef9c 100644 --- a/Source/PeleLMeX_Regrid.cpp +++ b/Source/PeleLMeX_Regrid.cpp @@ -316,8 +316,7 @@ PeleLM::MakeNewLevelFromCoarse( auto const dx = geom[lev].CellSizeArray(); Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); amrex::Print() << " with " << ba.numPts() << " cells, " << ba.size() - << " boxes," - << " over " + << " boxes," << " over " << static_cast(ba.numPts()) * vol / geom[0].ProbSize() * 100 << "% of the domain \n"; @@ -427,8 +426,7 @@ PeleLM::RemakeLevel( auto const dx = geom[lev].CellSizeArray(); Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); amrex::Print() << " with " << ba.numPts() << " cells," << ba.size() - << " boxes," - << " over " + << " boxes," << " over " << static_cast(ba.numPts()) * vol / geom[0].ProbSize() * 100 << "% of the domain \n"; diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index b7a508524..96673955b 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -121,6 +121,11 @@ PeleLM::Setup() amrex::Print() << " Using mixture-averaged transport with Soret effects" << std::endl; + if (m_soret_boundary_override) { + amrex::Print() << " Imposing inhomogenous Neumann conditions " + "for species on isothermal walls" + << std::endl; + } } } else { if (m_fixed_Le != 0) { @@ -361,6 +366,19 @@ PeleLM::readParameters() // diffusion ParmParse pptrans("transport"); pptrans.query("use_soret", m_use_soret); + 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; + } + } pp.query("use_wbar", m_use_wbar); pp.query("unity_Le", m_unity_Le); pp.query("fixed_Le", m_fixed_Le); From 145f8e05095b004dd8b84329ac22ad3f09cbbd85 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 21 Jun 2024 15:11:33 +0200 Subject: [PATCH 02/34] formatting --- Source/PeleLMeX_Diffusion.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 2927b10ba..a2ebd832d 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -176,8 +176,8 @@ PeleLM::computeDifferentialDiffusionTerms( #ifdef AMREX_USE_EB fluxDivergenceRD( GetVecOfConstPtrs(getSpeciesVect(a_time)), 0, GetVecOfPtrs(diffData->DT), - 0, GetVecOfArrOfPtrs(diffData->soret_fluxes), 0, {}, 0, NUM_SPECIES, intensiveFluxes, - 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), @@ -944,7 +944,7 @@ PeleLM::differentialDiffusionUpdate( const Box& ebx = mfi.nodaltilebox(idim); auto const& flux_soret = diffData->soret_fluxes[lev][idim].const_array(mfi); - auto const& flux_wbar = + auto const& flux_wbar = diffData->wbar_fluxes[lev][idim].const_array(mfi); auto const& rhoD_ec = beta_ec[idim].const_array(mfi); auto const& boundary_ar = spec_boundary[lev]->array(mfi); From ac393fd7926cf08d98eaccae07ae8a07586c9b0a Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 21 Jun 2024 15:14:11 +0200 Subject: [PATCH 03/34] formatting --- Source/PeleLMeX_Init.cpp | 3 ++- Source/PeleLMeX_Regrid.cpp | 6 ++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/Source/PeleLMeX_Init.cpp b/Source/PeleLMeX_Init.cpp index fbe0b14de..239b6c30c 100644 --- a/Source/PeleLMeX_Init.cpp +++ b/Source/PeleLMeX_Init.cpp @@ -38,7 +38,8 @@ PeleLM::MakeNewLevelFromScratch( auto const dx = geom[lev].CellSizeArray(); Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); amrex::Print() << " with " << ba.numPts() << " cells," << ba.size() - << " boxes," << " over " + << " boxes," + << " over " << static_cast(ba.numPts()) * vol / geom[0].ProbSize() * 100 << "% of the domain \n"; diff --git a/Source/PeleLMeX_Regrid.cpp b/Source/PeleLMeX_Regrid.cpp index b61a5ef9c..c4d3fb2ba 100644 --- a/Source/PeleLMeX_Regrid.cpp +++ b/Source/PeleLMeX_Regrid.cpp @@ -316,7 +316,8 @@ PeleLM::MakeNewLevelFromCoarse( auto const dx = geom[lev].CellSizeArray(); Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); amrex::Print() << " with " << ba.numPts() << " cells, " << ba.size() - << " boxes," << " over " + << " boxes," + << " over " << static_cast(ba.numPts()) * vol / geom[0].ProbSize() * 100 << "% of the domain \n"; @@ -426,7 +427,8 @@ PeleLM::RemakeLevel( auto const dx = geom[lev].CellSizeArray(); Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); amrex::Print() << " with " << ba.numPts() << " cells," << ba.size() - << " boxes," << " over " + << " boxes," + << " over " << static_cast(ba.numPts()) * vol / geom[0].ProbSize() * 100 << "% of the domain \n"; From de385ce525e1d5af992b42b0df830104cd78354c Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 21 Jun 2024 15:16:22 +0200 Subject: [PATCH 04/34] spelling --- Source/PeleLMeX_Setup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 96673955b..26b09631e 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -122,7 +122,7 @@ PeleLM::Setup() << " Using mixture-averaged transport with Soret effects" << std::endl; if (m_soret_boundary_override) { - amrex::Print() << " Imposing inhomogenous Neumann conditions " + amrex::Print() << " Imposing inhomogeneous Neumann conditions " "for species on isothermal walls" << std::endl; } From 174348cc84c724ed720a87ff20ce3e90b278e4b0 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 21 Jun 2024 17:51:28 +0200 Subject: [PATCH 05/34] add flag --- Source/PeleLMeX_Diffusion.cpp | 138 ++++++++++++++++---------------- Source/PeleLMeX_DiffusionOp.H | 2 +- Source/PeleLMeX_DiffusionOp.cpp | 5 +- 3 files changed, 73 insertions(+), 72 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index a2ebd832d..fa5a37677 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -306,7 +306,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( @@ -903,78 +902,79 @@ PeleLM::differentialDiffusionUpdate( // Get the species BCRec auto bcRecSpec = fetchBCRecArray(FIRSTSPEC, NUM_SPECIES); - Vector spec_boundary(finest_level + 1); - for (int lev = 0; lev <= finest_level; ++lev) { - auto* ldata_p = getLevelDataPtr(lev, AmrNewTime); - spec_boundary[lev] = new MultiFab( - ldata_p->state.boxArray(), ldata_p->state.DistributionMap(), NUM_SPECIES, - 1, MFInfo(), ldata_p->state.Factory()); - // copy the state into the boundaries for Dirichlet boundaries - MultiFab::Copy( - *spec_boundary[lev], ldata_p->state, FIRSTSPEC, 0, NUM_SPECIES, 1); - // to remain consistent with diffusion solve, we need to divide by density - // here if we have mix of Dirichlet and Isothermal conditions, we need to do - // work here, rather than letting nullptr do the work - for (int n = 0; n < NUM_SPECIES; n++) { - MultiFab::Divide(*spec_boundary[lev], ldata_p->state, DENSITY, n, 1, 1); - } - // no Isothermal conditions, just carry Dirichlet conditions through to - // diffusion solve - if (m_soret_boundary_override == 0) { - break; - } - // otherwise, lets overwrite the boundaries with the fluxes - 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); + Vector spec_boundary; + if (m_soret_boundary_override != 0) { + spec_boundary.resize(finest_level + 1); + 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); + } + + // 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(*spec_boundary[lev], 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 = - diffData->soret_fluxes[lev][idim].const_array(mfi); - auto const& flux_wbar = - diffData->wbar_fluxes[lev][idim].const_array(mfi); - auto const& rhoD_ec = beta_ec[idim].const_array(mfi); - auto const& boundary_ar = spec_boundary[lev]->array(mfi); - const int use_wbar = m_use_wbar; - amrex::ParallelFor( - ebx, - [flux_wbar, flux_soret, rhoD_ec, boundary_ar, idim, edomain, bc_lo, - bc_hi, use_wbar] 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) { - boundary_ar(idx[0], idx[1], idx[2], n) += - flux_wbar(i, j, k, n); + 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 = + diffData->soret_fluxes[lev][idim].const_array(mfi); + auto const& flux_wbar = + diffData->wbar_fluxes[lev][idim].const_array(mfi); + auto const& rhoD_ec = beta_ec[idim].const_array(mfi); + auto const& boundary_ar = 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] 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) { + 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); } - boundary_ar(idx[0], idx[1], idx[2], n) /= rhoD_ec(i, j, k, n); } - } - }); + }); + } } } } @@ -1022,7 +1022,7 @@ PeleLM::differentialDiffusionUpdate( GetVecOfConstPtrs( getDensityVect(AmrNewTime)), // this triggers proper scaling by density GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), 0, bcRecSpec, - NUM_SPECIES, 0, m_dt, spec_boundary); + NUM_SPECIES, 0, m_dt, GetVecOfConstPtrs(spec_boundary)); #endif // Add lagged Wbar term @@ -1261,7 +1261,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 diff --git a/Source/PeleLMeX_DiffusionOp.H b/Source/PeleLMeX_DiffusionOp.H index a365b01ab..465a16d82 100644 --- a/Source/PeleLMeX_DiffusionOp.H +++ b/Source/PeleLMeX_DiffusionOp.H @@ -37,7 +37,7 @@ public: int ncomp, int isPoissonSolve, amrex::Real dt, - amrex::Vector a_boundary = {nullptr}); + amrex::Vector const& a_boundary); #ifdef AMREX_USE_EB void diffuse_scalar( diff --git a/Source/PeleLMeX_DiffusionOp.cpp b/Source/PeleLMeX_DiffusionOp.cpp index 2309bbd92..78bf976a3 100644 --- a/Source/PeleLMeX_DiffusionOp.cpp +++ b/Source/PeleLMeX_DiffusionOp.cpp @@ -107,7 +107,7 @@ DiffusionOp::diffuse_scalar( int ncomp, int isPoissonSolve, Real a_dt, - Vector a_boundary) + Vector const& a_boundary) { BL_PROFILE("DiffusionOp::diffuse_scalar()"); @@ -117,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 @@ -231,7 +232,7 @@ 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); - if (a_boundary[lev] != nullptr) { + if (have_boundary != 0) { boundary.emplace_back( *a_boundary[lev], amrex::make_alias, comp, m_ncomp); } else { From f643ab2476f8044cf9a01701083188a081372a37 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 21 Jun 2024 17:58:26 +0200 Subject: [PATCH 06/34] fix for efield --- Source/PeleLMeX_Diffusion.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index fa5a37677..fec62f76c 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -993,7 +993,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,GetVecOfConstPtrs(spec_boundary)); // Ions one by one for (int n = 0; n < NUM_IONS; n++) { auto bcRecIons = fetchBCRecArray(FIRSTSPEC + NUM_SPECIES - NUM_IONS + n, 1); @@ -1006,7 +1006,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}} From 0cd9fab94fd5f6ce4d8a3142c8add627464159be Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 21 Jun 2024 17:59:57 +0200 Subject: [PATCH 07/34] formatting --- Source/PeleLMeX_Diffusion.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index fec62f76c..2ffeb64fd 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -993,7 +993,7 @@ PeleLM::differentialDiffusionUpdate( GetVecOfConstPtrs( getDensityVect(AmrNewTime)), // this triggers proper scaling by density GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), 0, bcRecSpec, - NUM_SPECIES - NUM_IONS, 0, m_dt,GetVecOfConstPtrs(spec_boundary)); + NUM_SPECIES - NUM_IONS, 0, m_dt, GetVecOfConstPtrs(spec_boundary)); // Ions one by one for (int n = 0; n < NUM_IONS; n++) { auto bcRecIons = fetchBCRecArray(FIRSTSPEC + NUM_SPECIES - NUM_IONS + n, 1); @@ -1006,7 +1006,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}} From c324e7271654e47d27801ffa380e467746198b61 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 21 Jun 2024 18:05:50 +0200 Subject: [PATCH 08/34] efield fix --- Source/Efield/PeleLMeX_EFPoisson.cpp | 2 +- Source/PeleLMeX_Diffusion.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Efield/PeleLMeX_EFPoisson.cpp b/Source/Efield/PeleLMeX_EFPoisson.cpp index f763f0daa..b18812f9e 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_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 2ffeb64fd..ed580a59d 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -1253,7 +1253,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 { From 22d13ff7b33323321fd4ecba320b71ebf2cd801f Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 21 Jun 2024 18:10:49 +0200 Subject: [PATCH 09/34] ignore eb for now --- Source/PeleLMeX_Diffusion.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index ed580a59d..2ffeb64fd 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -1253,7 +1253,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 { From 0f321fa5e6a9952c35022aea83b4a5590b2e709a Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 21 Jun 2024 18:11:47 +0200 Subject: [PATCH 10/34] formatting --- Source/Efield/PeleLMeX_EFPoisson.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Efield/PeleLMeX_EFPoisson.cpp b/Source/Efield/PeleLMeX_EFPoisson.cpp index b18812f9e..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, {}); } From 163dc645d32ca993daebb34a6bb4d6b31b5f16fc Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 21 Jun 2024 19:24:54 +0200 Subject: [PATCH 11/34] bool to int --- Source/PeleLMeX_Setup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 26b09631e..8dc863255 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -121,7 +121,7 @@ PeleLM::Setup() amrex::Print() << " Using mixture-averaged transport with Soret effects" << std::endl; - if (m_soret_boundary_override) { + if (m_soret_boundary_override != 0) { amrex::Print() << " Imposing inhomogeneous Neumann conditions " "for species on isothermal walls" << std::endl; From 22a8792c642e2068b70db6cfc786b432a529e5bc Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Mon, 24 Jun 2024 10:27:58 +0200 Subject: [PATCH 12/34] test case --- .../IsothermalSoretChannel/CMakeLists.txt | 8 ++ .../IsothermalSoretChannel/CMakeLists.txt~ | 8 ++ .../IsothermalSoretChannel/GNUmakefile | 37 ++++++ .../RegTests/IsothermalSoretChannel/README.md | 3 + .../RegTests/IsothermalSoretChannel/input.inp | 111 ++++++++++++++++ .../IsothermalSoretChannel/pelelmex_prob.H | 119 ++++++++++++++++++ .../IsothermalSoretChannel/pelelmex_prob.cpp | 16 +++ .../pelelmex_prob_parm.H | 18 +++ 8 files changed, 320 insertions(+) create mode 100644 Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt create mode 100644 Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt~ create mode 100644 Exec/RegTests/IsothermalSoretChannel/GNUmakefile create mode 100644 Exec/RegTests/IsothermalSoretChannel/README.md create mode 100644 Exec/RegTests/IsothermalSoretChannel/input.inp create mode 100644 Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H create mode 100644 Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp create mode 100644 Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H 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/CMakeLists.txt~ b/Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt~ new file mode 100644 index 000000000..cd74295e8 --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt~ @@ -0,0 +1,8 @@ +set(PELE_PHYSICS_EOS_MODEL Fuego) +set(PELE_PHYSICS_CHEMISTRY_MODEL drm19) +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..a4d3cb7f7 --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/input.inp @@ -0,0 +1,111 @@ +#---------------------- 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.Zst = 0.055 +peleLM.deltaT_crashIfFailing= 0 +#---------------------- 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 progress_variable + +#---------------------- AC CONTROL ------------------------------- +active_control.on = 0 # Use AC ? +# active_control.use_temp = 1 # Default in fuel mass, rather use iso-T position ? +# active_control.temperature = 1400.0 # Value of iso-T ? +# active_control.tau = 5.0e-4 # Control tau (should ~ 10 dt) +# active_control.height = 0.01 # Where is the flame held ? Default assumes coordinate along Y in 2D or Z in 3D. +# active_control.v = 1 # verbose +# active_control.method = 1 # Controller: 1 - Linear, 2 - Quadratic, 3 - Weighted quadratic +# active_control.velMax = 2.0 # Optional: limit inlet velocity +# active_control.changeMax = 0.1 # Optional: limit inlet velocity changes (absolute m/s) +# active_control.flow_dir = 1 # Optional: flame main direction. Default: AMREX_SPACEDIM-1 +# active_control.pseudo_gravity = 1 # Optional: add density proportional force to compensate for the acceleration + # of the gas due to inlet velocity changes + +#---------------------- 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 +peleLM.progressVariable.format = Cantera +peleLM.progressVariable.weights = H2O:1.0 +peleLM.progressVariable.coldState = H2O:0.0 +peleLM.progressVariable.hotState = H2O:0.012 + +#---------------------- 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 HR +amr.highT.max_level = 1 +amr.highT.value_greater = 500 +amr.highT.field_name = temp + +amr.gradT.max_level = 2 +amr.gradT.adjacent_difference_greater = 200 +amr.gradT.field_name = temp + +amr.HR.max_level = 2 +amr.HR.value_greater = 0.5e9 +amr.HR.field_name = HeatRelease + +#---------------------- 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..cb59fc174 --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H @@ -0,0 +1,119 @@ +#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_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}; + + // Define air + amrex::Real air_Y[NUM_SPECIES] = {0.0}; + air_Y[prob_parm.bathID] = 0.717; + air_Y[prob_parm.oxidID] = 0.183; + air_Y[prob_parm.fuelID] = 0.1; + + const amrex::Real T_low = 300.0; + const amrex::Real T_high = 700.0; + + state(i, j, k, TEMP) = (y - y_low) / (Ly) * (T_high - T_low) + T_low; + + for (int n = 0; n < NUM_SPECIES; ++n) { + massfrac[n] = air_Y[n]; + } + + state(i, j, k, VELX) = 10.0; + 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] = 300; + } else if (sgn == -1) { + s_ext[TEMP] = 700; + } + +} + +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..6c5eba670 --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp @@ -0,0 +1,16 @@ +#include +#include + +void +PeleLM::readProbParm() +{ + amrex::ParmParse pp("prob"); + + pp.query("P_mean", PeleLM::prob_parm->P_mean); + + // TODO: somewhat hard coded bath, fuel and oxid IDs + // should exist somewhere in PeleLM. + PeleLM::prob_parm->bathID = N2_ID; + PeleLM::prob_parm->fuelID = H2_ID; + PeleLM::prob_parm->oxidID = O2_ID; +} diff --git a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H new file mode 100644 index 000000000..e428ed3e1 --- /dev/null +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H @@ -0,0 +1,18 @@ +#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; + + int bathID{-1}; + int fuelID{-1}; + int oxidID{-1}; + +}; +#endif From 0b6613f47fc5ea74358a5fd3d231ad7660368272 Mon Sep 17 00:00:00 2001 From: Thomas Howarth <69794367+ThomasHowarth@users.noreply.github.com> Date: Mon, 24 Jun 2024 10:29:12 +0200 Subject: [PATCH 13/34] remove autosave file --- Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt~ | 8 -------- 1 file changed, 8 deletions(-) delete mode 100644 Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt~ diff --git a/Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt~ b/Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt~ deleted file mode 100644 index cd74295e8..000000000 --- a/Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt~ +++ /dev/null @@ -1,8 +0,0 @@ -set(PELE_PHYSICS_EOS_MODEL Fuego) -set(PELE_PHYSICS_CHEMISTRY_MODEL drm19) -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) From 7e47b23d24eac7cbfdc60813f792d1c9158cf4a1 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Mon, 24 Jun 2024 10:31:45 +0200 Subject: [PATCH 14/34] formatting --- Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H | 1 - Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H | 1 - 2 files changed, 2 deletions(-) diff --git a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H index cb59fc174..5c6b3743c 100644 --- a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H @@ -93,7 +93,6 @@ bcnormal( } else if (sgn == -1) { s_ext[TEMP] = 700; } - } AMREX_GPU_DEVICE diff --git a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H index e428ed3e1..009df563a 100644 --- a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H @@ -13,6 +13,5 @@ struct ProbParm int bathID{-1}; int fuelID{-1}; int oxidID{-1}; - }; #endif From 322a452cd9319a8dd9ef67aa4cd4638a9b388f89 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Mon, 24 Jun 2024 14:04:50 +0200 Subject: [PATCH 15/34] make test case more robust --- .../RegTests/IsothermalSoretChannel/input.inp | 44 +++++---------- .../IsothermalSoretChannel/pelelmex_prob.H | 53 ++++++++++++++----- .../IsothermalSoretChannel/pelelmex_prob.cpp | 24 +++++++-- .../pelelmex_prob_parm.H | 4 ++ 4 files changed, 76 insertions(+), 49 deletions(-) diff --git a/Exec/RegTests/IsothermalSoretChannel/input.inp b/Exec/RegTests/IsothermalSoretChannel/input.inp index a4d3cb7f7..d7360016c 100644 --- a/Exec/RegTests/IsothermalSoretChannel/input.inp +++ b/Exec/RegTests/IsothermalSoretChannel/input.inp @@ -22,8 +22,10 @@ amr.max_grid_size = 32 # max box size #---------------------- Problem ---------------------------------- prob.P_mean = 101325.0 -prob.Zst = 0.055 -peleLM.deltaT_crashIfFailing= 0 +prob.Vin = 10.0 +prob.Tlow = 400.0 +prob.Thigh = 800.0 +prob.phi = 0.5 #---------------------- Transport -------------------------------- transport.use_soret = 1 @@ -52,21 +54,7 @@ amr.dt_change_max = 1.1 # Maximum dt increase btw successive s #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 progress_variable - -#---------------------- AC CONTROL ------------------------------- -active_control.on = 0 # Use AC ? -# active_control.use_temp = 1 # Default in fuel mass, rather use iso-T position ? -# active_control.temperature = 1400.0 # Value of iso-T ? -# active_control.tau = 5.0e-4 # Control tau (should ~ 10 dt) -# active_control.height = 0.01 # Where is the flame held ? Default assumes coordinate along Y in 2D or Z in 3D. -# active_control.v = 1 # verbose -# active_control.method = 1 # Controller: 1 - Linear, 2 - Quadratic, 3 - Weighted quadratic -# active_control.velMax = 2.0 # Optional: limit inlet velocity -# active_control.changeMax = 0.1 # Optional: limit inlet velocity changes (absolute m/s) -# active_control.flow_dir = 1 # Optional: flame main direction. Default: AMREX_SPACEDIM-1 -# active_control.pseudo_gravity = 1 # Optional: add density proportional force to compensate for the acceleration - # of the gas due to inlet velocity changes +amr.derive_plot_vars = avg_pressure mag_vort mass_fractions mixture_fraction #---------------------- Derived CONTROLS ------------------------- peleLM.fuel_name = H2 @@ -74,37 +62,29 @@ peleLM.mixtureFraction.format = Cantera peleLM.mixtureFraction.type = mass peleLM.mixtureFraction.oxidTank = O2:0.233 N2:0.767 peleLM.mixtureFraction.fuelTank = H2:1.0 -peleLM.progressVariable.format = Cantera -peleLM.progressVariable.weights = H2O:1.0 -peleLM.progressVariable.coldState = H2O:0.0 -peleLM.progressVariable.hotState = H2O:0.012 #---------------------- 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. +#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 HR +amr.refinement_indicators = highT gradT amr.highT.max_level = 1 -amr.highT.value_greater = 500 +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 -amr.HR.max_level = 2 -amr.HR.value_greater = 0.5e9 -amr.HR.field_name = HeatRelease - #---------------------- Debug/HPC CONTROL------------------------- amrex.fpe_trap_invalid = 1 amrex.fpe_trap_zero = 1 diff --git a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H index 5c6b3743c..76ab5134d 100644 --- a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H @@ -10,6 +10,39 @@ #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; + amrex::Print() << prob_parm.fuelID << std::endl; + amrex::Print() << H2_ID << std::endl; + 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 @@ -38,22 +71,14 @@ pelelmex_initdata( auto eos = pele::physics::PhysicsType::eos(); amrex::Real massfrac[NUM_SPECIES] = {0.0}; - // Define air - amrex::Real air_Y[NUM_SPECIES] = {0.0}; - air_Y[prob_parm.bathID] = 0.717; - air_Y[prob_parm.oxidID] = 0.183; - air_Y[prob_parm.fuelID] = 0.1; + set_Y_from_Phi(prob_parm, massfrac); - const amrex::Real T_low = 300.0; - const amrex::Real T_high = 700.0; + 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; - for (int n = 0; n < NUM_SPECIES; ++n) { - massfrac[n] = air_Y[n]; - } - - state(i, j, k, VELX) = 10.0; + state(i, j, k, VELX) = prob_parm.Vin; state(i, j, k, VELY) = 0.0; amrex::Real rho_cgs, P_cgs; @@ -89,9 +114,9 @@ bcnormal( amrex::Real massfrac[NUM_SPECIES] = {0.0}; if (sgn == 1) { - s_ext[TEMP] = 300; + s_ext[TEMP] = prob_parm.Tlow; } else if (sgn == -1) { - s_ext[TEMP] = 700; + s_ext[TEMP] = prob_parm.Thigh; } } diff --git a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp index 6c5eba670..601464a2f 100644 --- a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp @@ -7,10 +7,28 @@ 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); - // TODO: somewhat hard coded bath, fuel and oxid IDs - // should exist somewhere in PeleLM. PeleLM::prob_parm->bathID = N2_ID; - PeleLM::prob_parm->fuelID = H2_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); + amrex::Print() << "H2_ID 1 " << H2_ID << std::endl; +#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 index 009df563a..6ae9189a8 100644 --- a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H @@ -9,6 +9,10 @@ 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}; From f6159d56ee06e539b0f11a8fd9049794d815492f Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Tue, 2 Jul 2024 18:30:33 +0200 Subject: [PATCH 16/34] remove prints --- Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H | 2 -- Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp | 1 - 2 files changed, 3 deletions(-) diff --git a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H index 76ab5134d..33c12a843 100644 --- a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H @@ -19,8 +19,6 @@ set_Y_from_Phi(ProbParm const& prob_parm, amrex::Real Y[]) auto eos = pele::physics::PhysicsType::eos(); amrex::Real Xt[NUM_SPECIES] = {0.0}; amrex::Real a = 0.0; - amrex::Print() << prob_parm.fuelID << std::endl; - amrex::Print() << H2_ID << std::endl; switch (prob_parm.fuelID) { #if defined(H2_ID) case H2_ID: diff --git a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp index 601464a2f..53c005d5c 100644 --- a/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp +++ b/Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp @@ -17,7 +17,6 @@ PeleLM::readProbParm() amrex::ParmParse pp_pele("peleLM"); std::string fuelName = ""; pp_pele.get("fuel_name", fuelName); - amrex::Print() << "H2_ID 1 " << H2_ID << std::endl; #if defined(H2_ID) if (fuelName == "H2") { PeleLM::prob_parm->fuelID = H2_ID; From d1869c7282880a83f7227dd9215be9ef68570a6f Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Tue, 2 Jul 2024 18:34:39 +0200 Subject: [PATCH 17/34] remove interaction with efield --- Source/PeleLMeX_Diffusion.cpp | 2 +- Source/PeleLMeX_Setup.cpp | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 2ffeb64fd..2992469b6 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -993,7 +993,7 @@ PeleLM::differentialDiffusionUpdate( GetVecOfConstPtrs( getDensityVect(AmrNewTime)), // this triggers proper scaling by density GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), 0, bcRecSpec, - NUM_SPECIES - NUM_IONS, 0, m_dt, GetVecOfConstPtrs(spec_boundary)); + 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); diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 8dc863255..00a7d60b6 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -377,6 +377,9 @@ PeleLM::readParameters() } if (isothermal) { m_soret_boundary_override = 1; +#if PELE_USE_EFIELD + amrex::Abort("Isothermal walls with Soret incompatible with Efield"); +#endif } } pp.query("use_wbar", m_use_wbar); From 7a92b9f507a27a7e27871ed103ab9473c0fc46d0 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 7 Aug 2024 13:32:34 +0200 Subject: [PATCH 18/34] fix for divu computation --- Source/PeleLMeX.H | 3 +- Source/PeleLMeX_Diffusion.cpp | 162 +++++++++++++++++++++++++++----- Source/PeleLMeX_DiffusionOp.H | 4 +- Source/PeleLMeX_DiffusionOp.cpp | 41 ++++++-- Source/PeleLMeX_K.H | 26 +++++ 5 files changed, 201 insertions(+), 35 deletions(-) diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index db69ffc4a..672d66b5a 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -583,7 +583,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 diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 2992469b6..2b22245ef 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( @@ -224,7 +226,89 @@ PeleLM::computeDifferentialDiffusionFluxes( // Species fluxes // 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); + int have_soret_fluxes = (a_soretfluxes.empty() || m_sdcIter == 0) ? 0 : 1; + int have_wbar_fluxes = (a_wbarfluxes.empty() || m_sdcIter == 0) ? 0 : 1; + for (int lev = 0; lev <= finest_level; ++lev) { + auto* ldata_p = getLevelDataPtr(lev, AmrOldTime); + spec_boundary[lev].define( + grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); + if (have_soret_fluxes == 0) { //no soret fluxes yet, set to zero and continue + spec_boundary[lev].setVal(0.0); + continue; + } + + // 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); + } + + // 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& rhoD_ec = beta_ec[idim].const_array(mfi); + auto const& flux_soret = (have_soret_fluxes != 0) ? a_soretfluxes[lev][idim]->const_array(mfi) : rhoD_ec; + auto const& flux_wbar = (have_wbar_fluxes != 0) ? a_wbarfluxes[lev][idim]->const_array(mfi) : rhoD_ec; + auto const& boundary_ar = spec_boundary[lev].array(mfi); + amrex::ParallelFor( + ebx, [flux_wbar, flux_soret, rhoD_ec, boundary_ar, idim, edomain, + bc_lo, bc_hi, have_wbar_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 (have_wbar_fluxes != 0) { + 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); + } + } + }); + } + } + } + } + + + #ifdef PELE_USE_EFIELD // Get the species diffusion fluxes from the DiffusionOp // Don't average down just yet @@ -234,7 +318,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); @@ -243,7 +327,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 @@ -254,7 +338,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 @@ -264,12 +348,12 @@ 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)); } } @@ -326,7 +410,7 @@ PeleLM::computeDifferentialDiffusionFluxes( 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 ) @@ -355,20 +439,27 @@ 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); + } 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 @@ -377,11 +468,32 @@ 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& gradWbar_arr = (have_boundary != 0) ? Wbar_boundary[lev].array(mfi) : Wbar_arr; + amrex::ParallelFor( - gbx, [rho_arr, rhoY_arr, - Wbar_arr] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - getMwmixGivenRY(i, j, k, rho_arr, rhoY_arr, Wbar_arr); - }); + gbx, [rho_arr, rhoY_arr, Wbar_arr, gradY_arr, gradWbar_arr, domain, have_boundary, phys_bc = m_phys_bc] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + getMwmixGivenRY(i, j, k, rho_arr, rhoY_arr, Wbar_arr); + if (have_boundary != 0) { //need to impose gradWbar on boundary for computeGradient + gradWbar_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,gradWbar_arr); + } + } + } + /* + std::cout << "Wbar" << std::endl; + std::cout << Wbar_arr(i,j,k) << std::endl; + std::cout << "gradWbar" << std::endl; + std::cout << gradWbar_arr(i,j,k) << std::endl; + */ + }); } } @@ -404,8 +516,8 @@ PeleLM::addWbarTerm( } } getDiffusionOp()->computeGradient( - GetVecOfArrOfPtrs(gradWbar), {}, // Don't need the laplacian out - GetVecOfConstPtrs(Wbar), bcRecSpec[0], do_avgDown); + GetVecOfArrOfPtrs(gradWbar), {}, // Don't need the laplacian out + GetVecOfConstPtrs(Wbar),GetVecOfConstPtrs(Wbar_boundary), bcRecSpec[0], do_avgDown); //------------------------------------------------------------------------ // add Wbar term to species fluxes @@ -529,8 +641,8 @@ PeleLM::addSoretTerm( } } getDiffusionOp()->computeGradient( - GetVecOfArrOfPtrs(gradT), {}, // Don't need the laplacian out - a_temp, bcRecTemp[0], do_avgDown); + GetVecOfArrOfPtrs(gradT), {}, // Don't need the laplacian out + a_temp, {}, bcRecTemp[0], do_avgDown); //------------------------------------------------------------------------ // add Soret term to species fluxes @@ -906,7 +1018,7 @@ PeleLM::differentialDiffusionUpdate( if (m_soret_boundary_override != 0) { spec_boundary.resize(finest_level + 1); for (int lev = 0; lev <= finest_level; ++lev) { - auto* ldata_p = getLevelDataPtr(lev, AmrNewTime); + auto* ldata_p = getLevelDataPtr(lev, AmrOldTime); spec_boundary[lev].define( grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); @@ -942,9 +1054,9 @@ PeleLM::differentialDiffusionUpdate( const Box& ebx = mfi.nodaltilebox(idim); auto const& flux_soret = diffData->soret_fluxes[lev][idim].const_array(mfi); - auto const& flux_wbar = - diffData->wbar_fluxes[lev][idim].const_array(mfi); auto const& rhoD_ec = beta_ec[idim].const_array(mfi); + auto const& flux_wbar = (m_use_wbar != 0) ? + diffData->wbar_fluxes[lev][idim].const_array(mfi) : rhoD_ec; auto const& boundary_ar = spec_boundary[lev].array(mfi); amrex::ParallelFor( ebx, [flux_wbar, flux_soret, rhoD_ec, boundary_ar, idim, edomain, @@ -1172,7 +1284,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 ) @@ -1415,7 +1527,7 @@ PeleLM::deltaTIter_update( 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 ) diff --git a/Source/PeleLMeX_DiffusionOp.H b/Source/PeleLMeX_DiffusionOp.H index 465a16d82..760122cd7 100644 --- a/Source/PeleLMeX_DiffusionOp.H +++ b/Source/PeleLMeX_DiffusionOp.H @@ -83,7 +83,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( @@ -108,6 +109,7 @@ public: 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 78bf976a3..26c0014b0 100644 --- a/Source/PeleLMeX_DiffusionOp.cpp +++ b/Source/PeleLMeX_DiffusionOp.cpp @@ -594,7 +594,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()"); @@ -609,6 +610,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 @@ -657,7 +659,8 @@ DiffusionOp::computeDiffFluxes( Vector> fluxes(finest_level + 1); Vector component; Vector laps; - + Vector boundary; + // Allow for component specific LinOp BC m_scal_apply_op->setDomainBC( m_pelelm->getDiffusionLinOpBC(Orientation::low, a_bcrec[comp]), @@ -669,6 +672,14 @@ 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 = { @@ -685,7 +696,8 @@ 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); @@ -838,6 +850,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 { @@ -845,7 +858,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); @@ -854,9 +866,10 @@ DiffusionOp::computeGradient( // Checks: one components only and 1 ghost cell at least AMREX_ASSERT(a_phi[0]->nComp() == 1); 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( m_pelelm->getDiffusionLinOpBC(Orientation::low, a_bcrec), @@ -865,13 +878,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()); - MultiFab::Copy(phi[lev], *a_phi[lev], 0, 0, 1, 1); - m_gradient_op->setLevelBC(lev, &phi[lev]); + 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); + + 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 20404845a..fba86ffe8 100644 --- a/Source/PeleLMeX_K.H +++ b/Source/PeleLMeX_K.H @@ -440,6 +440,32 @@ 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) noexcept +{ + using namespace amrex::literals; + + auto eos = pele::physics::PhysicsType::eos(); + 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 From e415fad3f285f061927c27196c61e2edd89a2132 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 7 Aug 2024 13:34:18 +0200 Subject: [PATCH 19/34] formatting --- Source/PeleLMeX.H | 2 +- Source/PeleLMeX_Diffusion.cpp | 153 ++++++++++++++++++-------------- Source/PeleLMeX_DiffusionOp.cpp | 12 ++- Source/PeleLMeX_K.H | 13 ++- 4 files changed, 100 insertions(+), 80 deletions(-) diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 672d66b5a..662f7e506 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -583,7 +583,7 @@ 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); /** diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 2b22245ef..b55fd303e 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -226,7 +226,7 @@ PeleLM::computeDifferentialDiffusionFluxes( // Species fluxes // 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); @@ -235,12 +235,13 @@ PeleLM::computeDifferentialDiffusionFluxes( for (int lev = 0; lev <= finest_level; ++lev) { auto* ldata_p = getLevelDataPtr(lev, AmrOldTime); spec_boundary[lev].define( - grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); - if (have_soret_fluxes == 0) { //no soret fluxes yet, set to zero and continue - spec_boundary[lev].setVal(0.0); - continue; + grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); + if (have_soret_fluxes == 0) { // no soret fluxes yet, set to zero and + // continue + spec_boundary[lev].setVal(0.0); + continue; } - + // 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 @@ -271,13 +272,20 @@ PeleLM::computeDifferentialDiffusionFluxes( const auto bc_hi = m_phys_bc.hi(idim); const Box& edomain = amrex::surroundingNodes(domain, idim); const Box& ebx = mfi.nodaltilebox(idim); - auto const& rhoD_ec = beta_ec[idim].const_array(mfi); - auto const& flux_soret = (have_soret_fluxes != 0) ? a_soretfluxes[lev][idim]->const_array(mfi) : rhoD_ec; - auto const& flux_wbar = (have_wbar_fluxes != 0) ? a_wbarfluxes[lev][idim]->const_array(mfi) : rhoD_ec; + auto const& rhoD_ec = beta_ec[idim].const_array(mfi); + auto const& flux_soret = + (have_soret_fluxes != 0) + ? a_soretfluxes[lev][idim]->const_array(mfi) + : rhoD_ec; + auto const& flux_wbar = (have_wbar_fluxes != 0) + ? a_wbarfluxes[lev][idim]->const_array(mfi) + : rhoD_ec; auto const& boundary_ar = spec_boundary[lev].array(mfi); amrex::ParallelFor( - ebx, [flux_wbar, flux_soret, rhoD_ec, boundary_ar, idim, edomain, - bc_lo, bc_hi, have_wbar_fluxes] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ebx, + [flux_wbar, flux_soret, rhoD_ec, boundary_ar, idim, edomain, bc_lo, + bc_hi, + have_wbar_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) && @@ -289,26 +297,24 @@ PeleLM::computeDifferentialDiffusionFluxes( 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 (have_wbar_fluxes != 0) { - 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); - } - } + 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 (have_wbar_fluxes != 0) { + 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); + } + } }); } } } } - - #ifdef PELE_USE_EFIELD // Get the species diffusion fluxes from the DiffusionOp // Don't average down just yet @@ -318,7 +324,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); @@ -327,7 +333,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 @@ -338,7 +344,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,GetVecOfConstPtrs(spec_boundary)); + do_avgDown, GetVecOfConstPtrs(spec_boundary)); #endif // Add the wbar term @@ -348,12 +354,14 @@ PeleLM::computeDifferentialDiffusionFluxes( addWbarTerm( a_fluxes, {}, GetVecOfConstPtrs(getSpeciesVect(a_time)), GetVecOfConstPtrs(getDensityVect(a_time)), - GetVecOfConstPtrs(getDiffusivityVect(a_time)),GetVecOfConstPtrs(spec_boundary)); + 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(spec_boundary)); + GetVecOfConstPtrs(getDiffusivityVect(a_time)), + GetVecOfConstPtrs(spec_boundary)); } } @@ -410,7 +418,7 @@ PeleLM::computeDifferentialDiffusionFluxes( 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 ) @@ -440,7 +448,7 @@ PeleLM::addWbarTerm( Vector const& a_spec, Vector const& a_rho, Vector const& a_beta, - Vector const& a_boundary) + Vector const& a_boundary) { //------------------------------------------------------------------------ // if a container for wbar fluxes is provided, fill it @@ -457,7 +465,8 @@ PeleLM::addWbarTerm( 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)); + Wbar_boundary[lev].define( + grids[lev], dmap[lev], 1, nGrow, MFInfo(), Factory(lev)); } const Box& domain = geom[lev].Domain(); #ifdef AMREX_USE_OMP @@ -468,32 +477,43 @@ 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& gradWbar_arr = (have_boundary != 0) ? Wbar_boundary[lev].array(mfi) : Wbar_arr; - + auto const& gradY_arr = + (have_boundary != 0) ? a_boundary[lev]->const_array(mfi) : Wbar_arr; + auto const& gradWbar_arr = + (have_boundary != 0) ? Wbar_boundary[lev].array(mfi) : Wbar_arr; + amrex::ParallelFor( - gbx, [rho_arr, rhoY_arr, Wbar_arr, gradY_arr, gradWbar_arr, domain, have_boundary, phys_bc = m_phys_bc] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - getMwmixGivenRY(i, j, k, rho_arr, rhoY_arr, Wbar_arr); - if (have_boundary != 0) { //need to impose gradWbar on boundary for computeGradient - gradWbar_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,gradWbar_arr); - } - } - } - /* - std::cout << "Wbar" << std::endl; - std::cout << Wbar_arr(i,j,k) << std::endl; - std::cout << "gradWbar" << std::endl; - std::cout << gradWbar_arr(i,j,k) << std::endl; - */ - }); + gbx, + [rho_arr, rhoY_arr, Wbar_arr, gradY_arr, gradWbar_arr, domain, + have_boundary, + phys_bc = m_phys_bc] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + getMwmixGivenRY(i, j, k, rho_arr, rhoY_arr, Wbar_arr); + if (have_boundary != 0) { // need to impose gradWbar on boundary for + // computeGradient + gradWbar_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, gradWbar_arr); + } + } + } + /* + std::cout << "Wbar" << std::endl; + std::cout << Wbar_arr(i,j,k) << std::endl; + std::cout << "gradWbar" << std::endl; + std::cout << gradWbar_arr(i,j,k) << std::endl; + */ + }); } } @@ -516,8 +536,9 @@ PeleLM::addWbarTerm( } } getDiffusionOp()->computeGradient( - GetVecOfArrOfPtrs(gradWbar), {}, // Don't need the laplacian out - GetVecOfConstPtrs(Wbar),GetVecOfConstPtrs(Wbar_boundary), bcRecSpec[0], do_avgDown); + GetVecOfArrOfPtrs(gradWbar), {}, // Don't need the laplacian out + GetVecOfConstPtrs(Wbar), GetVecOfConstPtrs(Wbar_boundary), bcRecSpec[0], + do_avgDown); //------------------------------------------------------------------------ // add Wbar term to species fluxes @@ -641,7 +662,7 @@ PeleLM::addSoretTerm( } } getDiffusionOp()->computeGradient( - GetVecOfArrOfPtrs(gradT), {}, // Don't need the laplacian out + GetVecOfArrOfPtrs(gradT), {}, // Don't need the laplacian out a_temp, {}, bcRecTemp[0], do_avgDown); //------------------------------------------------------------------------ @@ -1055,8 +1076,10 @@ PeleLM::differentialDiffusionUpdate( auto const& flux_soret = diffData->soret_fluxes[lev][idim].const_array(mfi); auto const& rhoD_ec = beta_ec[idim].const_array(mfi); - auto const& flux_wbar = (m_use_wbar != 0) ? - diffData->wbar_fluxes[lev][idim].const_array(mfi) : rhoD_ec; + auto const& flux_wbar = + (m_use_wbar != 0) + ? diffData->wbar_fluxes[lev][idim].const_array(mfi) + : rhoD_ec; auto const& boundary_ar = spec_boundary[lev].array(mfi); amrex::ParallelFor( ebx, [flux_wbar, flux_soret, rhoD_ec, boundary_ar, idim, edomain, @@ -1284,7 +1307,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 ) @@ -1527,7 +1550,7 @@ PeleLM::deltaTIter_update( 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 ) diff --git a/Source/PeleLMeX_DiffusionOp.cpp b/Source/PeleLMeX_DiffusionOp.cpp index 26c0014b0..e5fb0e12b 100644 --- a/Source/PeleLMeX_DiffusionOp.cpp +++ b/Source/PeleLMeX_DiffusionOp.cpp @@ -660,7 +660,7 @@ DiffusionOp::computeDiffFluxes( Vector component; Vector laps; Vector boundary; - + // Allow for component specific LinOp BC m_scal_apply_op->setDomainBC( m_pelelm->getDiffusionLinOpBC(Orientation::low, a_bcrec[comp]), @@ -679,7 +679,6 @@ DiffusionOp::computeDiffFluxes( boundary.emplace_back(phi[lev], amrex::make_alias, comp, m_ncomp); } - int doZeroVisc = 1; int addTurbContrib = 1; Vector subBCRec = { @@ -697,7 +696,6 @@ DiffusionOp::computeDiffFluxes( m_scal_apply_op->setBCoeffs(lev, GetArrOfConstPtrs(bcoeff_ec)); #endif m_scal_apply_op->setLevelBC(lev, &boundary[lev]); - } MLMG mlmg(*m_scal_apply_op); @@ -866,10 +864,10 @@ DiffusionOp::computeGradient( // Checks: one components only and 1 ghost cell at least AMREX_ASSERT(a_phi[0]->nComp() == 1); 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( m_pelelm->getDiffusionLinOpBC(Orientation::low, a_bcrec), @@ -888,14 +886,14 @@ DiffusionOp::computeGradient( 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); + MultiFab::Copy(phi[lev], *a_phi[lev], 0, 0, 1, 1); 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); diff --git a/Source/PeleLMeX_K.H b/Source/PeleLMeX_K.H index fba86ffe8..e683721ff 100644 --- a/Source/PeleLMeX_K.H +++ b/Source/PeleLMeX_K.H @@ -448,7 +448,7 @@ getGradMwmixGivengradYMwmix( int j, int k, amrex::Array4 const& gradY, - amrex::Array4 const& Mwmix, + amrex::Array4 const& Mwmix, amrex::Array4 const& gradMwmix) noexcept { using namespace amrex::literals; @@ -456,16 +456,15 @@ getGradMwmixGivengradYMwmix( auto eos = pele::physics::PhysicsType::eos(); amrex::Real imw[NUM_SPECIES] = {0.0_rt}; eos.inv_molecular_weight(imw); - gradMwmix(i,j,k) = 0.0; + 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]; + 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); + // 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 From 44f85e3bbe532f06adaccda6a6ca82dcaf4b201d Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 7 Aug 2024 13:36:34 +0200 Subject: [PATCH 20/34] remove prints --- Source/PeleLMeX_Diffusion.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index b55fd303e..cfadd7700 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -507,12 +507,6 @@ PeleLM::addWbarTerm( } } } - /* - std::cout << "Wbar" << std::endl; - std::cout << Wbar_arr(i,j,k) << std::endl; - std::cout << "gradWbar" << std::endl; - std::cout << gradWbar_arr(i,j,k) << std::endl; - */ }); } } From db1814bb60e56cbe7ef7168c22ce7c79dc718a62 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 7 Aug 2024 13:54:48 +0200 Subject: [PATCH 21/34] EB compatibility and Efield fix --- Source/Efield/PeleLMeX_EFIonDrift.cpp | 4 ++-- Source/Efield/PeleLMeX_EFNLSolve.cpp | 4 ++-- Source/PeleLMeX_Diffusion.cpp | 2 +- Source/PeleLMeX_DiffusionOp.cpp | 29 +++++++++++++++++++++------ 4 files changed, 28 insertions(+), 11 deletions(-) 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/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index cfadd7700..febeb0655 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -1382,7 +1382,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 { diff --git a/Source/PeleLMeX_DiffusionOp.cpp b/Source/PeleLMeX_DiffusionOp.cpp index e5fb0e12b..7be921e9e 100644 --- a/Source/PeleLMeX_DiffusionOp.cpp +++ b/Source/PeleLMeX_DiffusionOp.cpp @@ -326,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()"); @@ -336,7 +337,8 @@ 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 AMREX_ASSERT(m_ncomp == 1 || m_ncomp == ncomp); @@ -413,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( @@ -444,7 +447,12 @@ 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]); } @@ -735,7 +743,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()"); @@ -750,6 +759,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 @@ -800,7 +810,8 @@ DiffusionOp::computeDiffFluxes( Vector> ebfluxes; Vector component; Vector laps; - + Vector boundary; + // Allow for component specific LinOp BC m_scal_apply_op->setDomainBC( m_pelelm->getDiffusionLinOpBC(Orientation::low, a_bcrec[comp]), @@ -814,6 +825,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}; @@ -824,7 +841,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]); } From 236daea1c1027fd9bb0c9ad91da322094624aaae Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 7 Aug 2024 13:57:02 +0200 Subject: [PATCH 22/34] formatting --- Source/PeleLMeX_DiffusionOp.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Source/PeleLMeX_DiffusionOp.cpp b/Source/PeleLMeX_DiffusionOp.cpp index 7be921e9e..9ac11506d 100644 --- a/Source/PeleLMeX_DiffusionOp.cpp +++ b/Source/PeleLMeX_DiffusionOp.cpp @@ -338,7 +338,7 @@ DiffusionOp::diffuse_scalar( 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 AMREX_ASSERT(m_ncomp == 1 || m_ncomp == ncomp); @@ -448,7 +448,8 @@ DiffusionOp::diffuse_scalar( rhs.emplace_back( *a_rhs[lev], amrex::make_alias, rhs_comp + comp, m_ncomp); if (have_boundary != 0) { - boundary.emplace_back(*a_boundary[lev], amrex::make_alias, comp, m_ncomp); + boundary.emplace_back( + *a_boundary[lev], amrex::make_alias, comp, m_ncomp); } else { boundary.emplace_back(phi[lev], amrex::make_alias, comp, m_ncomp); } @@ -811,7 +812,7 @@ DiffusionOp::computeDiffFluxes( Vector component; Vector laps; Vector boundary; - + // Allow for component specific LinOp BC m_scal_apply_op->setDomainBC( m_pelelm->getDiffusionLinOpBC(Orientation::low, a_bcrec[comp]), From 478467ef152b211185c321e174158613c712e6bd Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 7 Aug 2024 15:29:57 +0200 Subject: [PATCH 23/34] header fix --- Source/PeleLMeX_DiffusionOp.H | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Source/PeleLMeX_DiffusionOp.H b/Source/PeleLMeX_DiffusionOp.H index 760122cd7..d4705b37d 100644 --- a/Source/PeleLMeX_DiffusionOp.H +++ b/Source/PeleLMeX_DiffusionOp.H @@ -59,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( @@ -102,7 +103,8 @@ 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( From d439b9359e84f16126f25acf6a51eb5fa274a3b6 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 7 Aug 2024 15:39:24 +0200 Subject: [PATCH 24/34] missing argument --- Source/PeleLMeX_Diffusion.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index febeb0655..046c43bed 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -411,7 +411,7 @@ 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 { @@ -1293,7 +1293,7 @@ PeleLM::differentialDiffusionUpdate( GetVecOfConstPtrs(getTempVect(AmrNewTime)), 0, {}, GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), NUM_SPECIES, GetVecOfConstPtrs(EBvalue), GetVecOfConstPtrs(EBdiff), bcRecTemp, 1, - do_avgDown); + do_avgDown, {}); } else #endif { @@ -1537,7 +1537,7 @@ PeleLM::deltaTIter_update( GetVecOfConstPtrs(getTempVect(AmrNewTime)), 0, {}, GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), NUM_SPECIES, GetVecOfConstPtrs(EBvalue), GetVecOfConstPtrs(EBdiff), bcRecTemp, 1, - do_avgDown); + do_avgDown, {}); } else #endif { From 3c7d3edeff24ed9403d2612fbb6e019861b1027f Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 18 Sep 2024 14:49:54 +0200 Subject: [PATCH 25/34] explicit flux calculation, plus wbar fix --- Source/PeleLMeX.H | 8 + Source/PeleLMeX_Diffusion.cpp | 366 +++++++++++++++++++--------------- 2 files changed, 209 insertions(+), 165 deletions(-) diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index fdaad1918..256eae737 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -552,6 +552,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, diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 046c43bed..dedb64d76 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -206,6 +206,171 @@ PeleLM::computeDifferentialDiffusionTerms( #endif } +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(getSpeciesVect(a_time)), + GetVecOfConstPtrs(getTempVect(a_time)), + GetVecOfConstPtrs(getDiffusivityVect(a_time))); + } else { // have the lagged ones, alias to 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); + } + } + }); + } + } + } + if (need_explicit_fluxes && m_use_wbar != 0) { // need to iterate to get the + // proper wbar flux + Vector> wbarfluxes(finest_level + 1); + for (int lev = 0; lev <= finest_level; lev++) { + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { + wbarfluxes[lev][idim].define( + grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); + wbarfluxes[lev][idim].setVal(0.0); + } + } + int iter_max = 10; + for (int iter = 0; iter < iter_max; iter++) { + // based on existing flux, compute required wbar flux + addWbarTerm( + GetVecOfArrOfPtrs(wbarfluxes), GetVecOfArrOfPtrs(wbarfluxes), + GetVecOfConstPtrs(getSpeciesVect(a_time)), + GetVecOfConstPtrs(getDensityVect(a_time)), + GetVecOfConstPtrs(getDiffusivityVect(a_time)), + GetVecOfConstPtrs(a_spec_boundary)); + for (int lev = 0; lev <= finest_level; ++lev) { + auto* ldata_p = getLevelDataPtr(lev, a_time); + + 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& rhoD_ec = beta_ec[idim].const_array(mfi); + auto const& flux_wbar = wbarfluxes[lev][idim].const_array(mfi); + auto const& flux_soret = soretfluxes[lev][idim]->const_array(mfi); + auto const& boundary_ar = a_spec_boundary[lev]->array(mfi); + amrex::ParallelFor( + ebx, [flux_soret, flux_wbar, rhoD_ec, boundary_ar, idim, edomain, + bc_lo, bc_hi, + iter] 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) + flux_wbar(i, j, k, n)) / + rhoD_ec(i, j, k, n); + } + } + }); + } + } + } + } + } +} + void PeleLM::computeDifferentialDiffusionFluxes( const TimeStamp& a_time, @@ -223,24 +388,15 @@ PeleLM::computeDifferentialDiffusionFluxes( BL_PROFILE("PeleLMeX::computeDifferentialDiffusionFluxes()"); //---------------------------------------------------------------- - // Species fluxes - // 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); - int have_soret_fluxes = (a_soretfluxes.empty() || m_sdcIter == 0) ? 0 : 1; - int have_wbar_fluxes = (a_wbarfluxes.empty() || m_sdcIter == 0) ? 0 : 1; + // this is the same regardless of lagged or not for (int lev = 0; lev <= finest_level; ++lev) { - auto* ldata_p = getLevelDataPtr(lev, AmrOldTime); + auto* ldata_p = getLevelDataPtr(lev, a_time); spec_boundary[lev].define( grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); - if (have_soret_fluxes == 0) { // no soret fluxes yet, set to zero and - // continue - spec_boundary[lev].setVal(0.0); - continue; - } // 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 @@ -251,69 +407,14 @@ PeleLM::computeDifferentialDiffusionFluxes( for (int n = 0; n < NUM_SPECIES; n++) { MultiFab::Divide(spec_boundary[lev], ldata_p->state, DENSITY, n, 1, 1); } - - // 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& rhoD_ec = beta_ec[idim].const_array(mfi); - auto const& flux_soret = - (have_soret_fluxes != 0) - ? a_soretfluxes[lev][idim]->const_array(mfi) - : rhoD_ec; - auto const& flux_wbar = (have_wbar_fluxes != 0) - ? a_wbarfluxes[lev][idim]->const_array(mfi) - : rhoD_ec; - auto const& boundary_ar = spec_boundary[lev].array(mfi); - amrex::ParallelFor( - ebx, - [flux_wbar, flux_soret, rhoD_ec, boundary_ar, idim, edomain, bc_lo, - bc_hi, - have_wbar_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 (have_wbar_fluxes != 0) { - 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); - } - } - }); - } - } } + // 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); #ifdef PELE_USE_EFIELD // Get the species diffusion fluxes from the DiffusionOp @@ -366,6 +467,7 @@ PeleLM::computeDifferentialDiffusionFluxes( } // Add the Soret term + // std::cout << "Calculating afterwards" << std::endl; if (m_use_soret != 0) { int need_soret_fluxes = (a_soretfluxes.empty()) ? 0 : 1; if (need_soret_fluxes == 0) { @@ -604,17 +706,24 @@ PeleLM::addWbarTerm( for (int n = 0; n < NUM_SPECIES; n++) { y[n] = rhoY(i, j, k, n) * rho_inv; } - amrex::Real WBAR = 0.0; - eos.Y2WBAR(y, WBAR); - WBAR *= 0.001; + amrex::Real imw[NUM_SPECIES] = {0.0}; + eos.inv_molecular_weight(imw); + // amrex::Real WBAR = 0.0; + // eos.Y2WBAR(y, WBAR); + // WBAR *= 0.001; for (int n = 0; n < NUM_SPECIES; n++) { + // imw[n] *= 1000.0; + // spFlux_ar(i, j, k, n) -= + // y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); spFlux_ar(i, j, k, n) -= - y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); + y[n] * imw[n] * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); } if (need_wbar_fluxes != 0) { for (int n = 0; n < NUM_SPECIES; n++) { + // spwbarFlux_ar(i, j, k, n) = + // -y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); spwbarFlux_ar(i, j, k, n) = - -y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); + -y[n] * imw[n] * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); } } }); @@ -675,36 +784,14 @@ PeleLM::addSoretTerm( #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif { - FArrayBox rhoY_ed; + // 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); @@ -739,12 +826,14 @@ 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 + Real dev = 0.0; amrex::ParallelFor( - ebx, - [need_soret_fluxes, gradT_ar, beta_ar, T, spFlux_ar, - spsoretFlux_ar] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ebx, [need_soret_fluxes, gradT_ar, beta_ar, T, spFlux_ar, domain, + spsoretFlux_ar, phys_bc = m_phys_bc, + dev] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int idx[3] = {i, j, k}; for (int n = 0; n < NUM_SPECIES; n++) { spFlux_ar(i, j, k, n) -= beta_ar(i, j, k, n) * gradT_ar(i, j, k) / T(i, j, k); @@ -1032,8 +1121,9 @@ PeleLM::differentialDiffusionUpdate( 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, AmrOldTime); + auto* ldata_p = getLevelDataPtr(lev, AmrNewTime); spec_boundary[lev].define( grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); @@ -1046,66 +1136,12 @@ PeleLM::differentialDiffusionUpdate( for (int n = 0; n < NUM_SPECIES; n++) { MultiFab::Divide(spec_boundary[lev], ldata_p->state, DENSITY, n, 1, 1); } - - // 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 = - diffData->soret_fluxes[lev][idim].const_array(mfi); - auto const& rhoD_ec = beta_ec[idim].const_array(mfi); - auto const& flux_wbar = - (m_use_wbar != 0) - ? diffData->wbar_fluxes[lev][idim].const_array(mfi) - : rhoD_ec; - auto const& boundary_ar = 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] 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) { - 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); - } - } - }); - } - } } + // 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 From cea2e6309f95a4e8cfa66bb13bb5dad99173cb8f Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 18 Sep 2024 14:57:14 +0200 Subject: [PATCH 26/34] removed unnecessary species call in soret term --- Source/PeleLMeX.H | 1 - Source/PeleLMeX_Diffusion.cpp | 11 +++-------- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 256eae737..9b16d9fc6 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -606,7 +606,6 @@ public: a_spfluxes, const amrex::Vector>& a_spsoretfluxes, - amrex::Vector const& a_spec, amrex::Vector const& a_temp, amrex::Vector const& a_beta); diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index dedb64d76..0920b82ca 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -227,8 +227,7 @@ PeleLM::correctIsothermalBoundary( } } addSoretTerm( - soretfluxes, soretfluxes, GetVecOfConstPtrs(getSpeciesVect(a_time)), - GetVecOfConstPtrs(getTempVect(a_time)), + soretfluxes, soretfluxes, GetVecOfConstPtrs(getTempVect(a_time)), GetVecOfConstPtrs(getDiffusivityVect(a_time))); } else { // have the lagged ones, alias to to them for (int lev = 0; lev <= finest_level; lev++) { @@ -472,13 +471,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))); } } @@ -737,7 +734,6 @@ void PeleLM::addSoretTerm( const Vector>& a_spfluxes, const Vector>& a_spsoretfluxes, - Vector const& a_spec, Vector const& a_temp, Vector const& a_beta) { @@ -833,7 +829,6 @@ PeleLM::addSoretTerm( ebx, [need_soret_fluxes, gradT_ar, beta_ar, T, spFlux_ar, domain, spsoretFlux_ar, phys_bc = m_phys_bc, dev] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - int idx[3] = {i, j, k}; for (int n = 0; n < NUM_SPECIES; n++) { spFlux_ar(i, j, k, n) -= beta_ar(i, j, k, n) * gradT_ar(i, j, k) / T(i, j, k); From fd18c3285df284702f799c2f0796a2558ae65eee Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 18 Sep 2024 15:09:10 +0200 Subject: [PATCH 27/34] remove unnecessary captures --- Source/PeleLMeX_Diffusion.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 0920b82ca..d7660c660 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -340,9 +340,9 @@ PeleLM::correctIsothermalBoundary( auto const& flux_soret = soretfluxes[lev][idim]->const_array(mfi); auto const& boundary_ar = a_spec_boundary[lev]->array(mfi); amrex::ParallelFor( - ebx, [flux_soret, flux_wbar, rhoD_ec, boundary_ar, idim, edomain, - bc_lo, bc_hi, - iter] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ebx, + [flux_soret, flux_wbar, rhoD_ec, boundary_ar, idim, edomain, + bc_lo, bc_hi] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { int idx[3] = {i, j, k}; bool on_lo = (bc_lo == BoundaryCondition::BCNoSlipWallIsotherm || @@ -826,9 +826,9 @@ PeleLM::addSoretTerm( // with beta_m = rho * D_m * chi_m below Real dev = 0.0; amrex::ParallelFor( - ebx, [need_soret_fluxes, gradT_ar, beta_ar, T, spFlux_ar, domain, - spsoretFlux_ar, phys_bc = m_phys_bc, - dev] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ebx, + [need_soret_fluxes, gradT_ar, beta_ar, T, spFlux_ar, + spsoretFlux_ar] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { for (int n = 0; n < NUM_SPECIES; n++) { spFlux_ar(i, j, k, n) -= beta_ar(i, j, k, n) * gradT_ar(i, j, k) / T(i, j, k); From 3178d306d6a4b0882faa6d36c20028cf999cf19d Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 18 Sep 2024 15:14:54 +0200 Subject: [PATCH 28/34] unnecessary variable --- Source/PeleLMeX_Diffusion.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index d7660c660..cfeb82888 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -824,7 +824,6 @@ PeleLM::addSoretTerm( // Soret flux is : - rho * D_m * chi_m * \nabla T / T // with beta_m = rho * D_m * chi_m below - Real dev = 0.0; amrex::ParallelFor( ebx, [need_soret_fluxes, gradT_ar, beta_ar, T, spFlux_ar, From 30c0f08f2f435cdf33252092a84d59c84b5689a3 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 18 Sep 2024 15:23:36 +0200 Subject: [PATCH 29/34] remove commented print --- Source/PeleLMeX_Diffusion.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index cfeb82888..c2fc3880e 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -466,7 +466,6 @@ PeleLM::computeDifferentialDiffusionFluxes( } // Add the Soret term - // std::cout << "Calculating afterwards" << std::endl; if (m_use_soret != 0) { int need_soret_fluxes = (a_soretfluxes.empty()) ? 0 : 1; if (need_soret_fluxes == 0) { From 5011524f1182c041e264ae392b683d251ac4b3f5 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 18 Sep 2024 16:04:59 +0200 Subject: [PATCH 30/34] go back to wbar for wbarflux --- Source/PeleLMeX_Diffusion.cpp | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index c2fc3880e..e6e685095 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -702,24 +702,17 @@ PeleLM::addWbarTerm( for (int n = 0; n < NUM_SPECIES; n++) { y[n] = rhoY(i, j, k, n) * rho_inv; } - amrex::Real imw[NUM_SPECIES] = {0.0}; - eos.inv_molecular_weight(imw); - // amrex::Real WBAR = 0.0; - // eos.Y2WBAR(y, WBAR); - // WBAR *= 0.001; + amrex::Real WBAR = 0.0; + eos.Y2WBAR(y, WBAR); + WBAR *= 0.001; for (int n = 0; n < NUM_SPECIES; n++) { - // imw[n] *= 1000.0; - // spFlux_ar(i, j, k, n) -= - // y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); - spFlux_ar(i, j, k, n) -= - y[n] * imw[n] * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); + spFlux_ar(i, j, k, n) -= + y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); } if (need_wbar_fluxes != 0) { for (int n = 0; n < NUM_SPECIES; n++) { - // spwbarFlux_ar(i, j, k, n) = - // -y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); - spwbarFlux_ar(i, j, k, n) = - -y[n] * imw[n] * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); + spwbarFlux_ar(i, j, k, n) = + -y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); } } }); From dd5699dbd0bd4c67c09d0c8117d6f6b88d9adcf5 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Wed, 18 Sep 2024 17:37:05 +0200 Subject: [PATCH 31/34] formatting --- Source/PeleLMeX_Diffusion.cpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index e6e685095..b45261b94 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -308,7 +308,7 @@ PeleLM::correctIsothermalBoundary( } int iter_max = 10; for (int iter = 0; iter < iter_max; iter++) { - // based on existing flux, compute required wbar flux + // based on existing gradient, compute the wbar flux addWbarTerm( GetVecOfArrOfPtrs(wbarfluxes), GetVecOfArrOfPtrs(wbarfluxes), GetVecOfConstPtrs(getSpeciesVect(a_time)), @@ -317,7 +317,6 @@ PeleLM::correctIsothermalBoundary( GetVecOfConstPtrs(a_spec_boundary)); for (int lev = 0; lev <= finest_level; ++lev) { auto* ldata_p = getLevelDataPtr(lev, a_time); - MultiFab& ldata_beta_cc = ldata_p->diff_cc; const Box& domain = geom[lev].Domain(); int doZeroVisc = 1; @@ -685,8 +684,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, @@ -702,17 +702,17 @@ PeleLM::addWbarTerm( for (int n = 0; n < NUM_SPECIES; n++) { y[n] = rhoY(i, j, k, n) * rho_inv; } - amrex::Real WBAR = 0.0; - eos.Y2WBAR(y, WBAR); - WBAR *= 0.001; + amrex::Real WBAR = 0.0; + eos.Y2WBAR(y, WBAR); + WBAR *= 0.001; for (int n = 0; n < NUM_SPECIES; n++) { - spFlux_ar(i, j, k, n) -= - y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); + spFlux_ar(i, j, k, n) -= + y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); } if (need_wbar_fluxes != 0) { for (int n = 0; n < NUM_SPECIES; n++) { - spwbarFlux_ar(i, j, k, n) = - -y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); + spwbarFlux_ar(i, j, k, n) = + -y[n] / WBAR * beta_ar(i, j, k, n) * gradWbar_ar(i, j, k); } } }); @@ -783,7 +783,7 @@ PeleLM::addSoretTerm( // 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); From da1d363eec5cc0bcbad5c92e049a64a12518beb3 Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Thu, 19 Sep 2024 18:15:34 +0200 Subject: [PATCH 32/34] CGS tweak and residual calculation --- Source/PeleLMeX_Diffusion.cpp | 42 ++++++++++++++++++++++++++--------- Source/PeleLMeX_K.H | 2 +- Source/PeleLMeX_Setup.cpp | 2 +- 3 files changed, 33 insertions(+), 13 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index b45261b94..0736819eb 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -217,6 +217,7 @@ PeleLM::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++) { @@ -237,12 +238,10 @@ PeleLM::correctIsothermalBoundary( } } } - 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(); @@ -296,17 +295,24 @@ PeleLM::correctIsothermalBoundary( } } } + if (need_explicit_fluxes && m_use_wbar != 0) { // need to iterate to get the // proper wbar flux + // check that we really don't have any wbar fluxes coming in + AMREX_ALWAYS_ASSERT(a_wbarfluxes.empty()); Vector> wbarfluxes(finest_level + 1); + Vector residual(finest_level + 1); + const int iter_max = 10; + const Real tol = 1e-10; for (int lev = 0; lev <= finest_level; lev++) { + // residual over all species + residual[lev].define(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { wbarfluxes[lev][idim].define( grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); wbarfluxes[lev][idim].setVal(0.0); } } - int iter_max = 10; for (int iter = 0; iter < iter_max; iter++) { // based on existing gradient, compute the wbar flux addWbarTerm( @@ -316,6 +322,7 @@ PeleLM::correctIsothermalBoundary( GetVecOfConstPtrs(getDiffusivityVect(a_time)), GetVecOfConstPtrs(a_spec_boundary)); for (int lev = 0; lev <= finest_level; ++lev) { + residual[lev].setVal(0.0); auto* ldata_p = getLevelDataPtr(lev, a_time); MultiFab& ldata_beta_cc = ldata_p->diff_cc; const Box& domain = geom[lev].Domain(); @@ -327,7 +334,7 @@ PeleLM::correctIsothermalBoundary( #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(ldata_beta_cc, TilingIfNotGPU()); mfi.isValid(); + for (MFIter mfi(*a_spec_boundary[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { const auto bc_lo = m_phys_bc.lo(idim); @@ -338,10 +345,12 @@ PeleLM::correctIsothermalBoundary( auto const& flux_wbar = wbarfluxes[lev][idim].const_array(mfi); auto const& flux_soret = soretfluxes[lev][idim]->const_array(mfi); auto const& boundary_ar = a_spec_boundary[lev]->array(mfi); + auto const& residual_ar = residual[lev].array(mfi); amrex::ParallelFor( ebx, [flux_soret, flux_wbar, rhoD_ec, boundary_ar, idim, edomain, - bc_lo, bc_hi] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + bc_lo, bc_hi, + residual_ar] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { int idx[3] = {i, j, k}; bool on_lo = (bc_lo == BoundaryCondition::BCNoSlipWallIsotherm || @@ -356,15 +365,26 @@ PeleLM::correctIsothermalBoundary( idx[idim] -= 1; } for (int n = 0; n < NUM_SPECIES; n++) { + // hold the old value + Real b_old = boundary_ar(idx[0], idx[1], idx[2], n); + // update the value boundary_ar(idx[0], idx[1], idx[2], n) = (flux_soret(i, j, k, n) + flux_wbar(i, j, k, n)) / rhoD_ec(i, j, k, n); + // compute residual error + Real res = + std::abs(boundary_ar(idx[0], idx[1], idx[2], n) - b_old); + residual_ar(idx[0], idx[1], idx[2]) = + std::max(residual_ar(idx[0], idx[1], idx[2]), res); } } }); } } } + if (residual[finest_level].norm0(0, 1) < tol) { + break; + } } } } @@ -576,18 +596,19 @@ PeleLM::addWbarTerm( 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& gradWbar_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, gradY_arr, gradWbar_arr, domain, + [rho_arr, rhoY_arr, Wbar_arr, gradY_arr, Wbar_boundary_arr, domain, have_boundary, phys_bc = m_phys_bc] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { getMwmixGivenRY(i, j, k, rho_arr, rhoY_arr, Wbar_arr); if (have_boundary != 0) { // need to impose gradWbar on boundary for - // computeGradient - gradWbar_arr(i, j, k) = Wbar_arr(i, j, k); + // computeGradient + // for dirichlet boundaries, we'll overwrite inhomog neumann ones + 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); @@ -600,7 +621,7 @@ PeleLM::addWbarTerm( (idx[idim] > domain.bigEnd(idim)); if (on_lo || on_hi) { getGradMwmixGivengradYMwmix( - i, j, k, gradY_arr, Wbar_arr, gradWbar_arr); + i, j, k, gradY_arr, Wbar_arr, Wbar_boundary_arr); } } } @@ -772,7 +793,6 @@ 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++) { diff --git a/Source/PeleLMeX_K.H b/Source/PeleLMeX_K.H index e683721ff..4c234ca7a 100644 --- a/Source/PeleLMeX_K.H +++ b/Source/PeleLMeX_K.H @@ -458,7 +458,7 @@ getGradMwmixGivengradYMwmix( 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 + // 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) diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 88056bc50..81d8d153f 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -406,6 +406,7 @@ PeleLM::readParameters() // diffusion 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++) { @@ -422,7 +423,6 @@ PeleLM::readParameters() #endif } } - pp.query("use_wbar", m_use_wbar); pp.query("unity_Le", m_unity_Le); pp.query("fixed_Le", m_fixed_Le); pp.query("fixed_Pr", m_fixed_Pr); From 1c094cad18d0963c9fa483d594d0de72d6b44a9b Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 15 Nov 2024 13:57:25 +0100 Subject: [PATCH 33/34] conflict fix --- Source/PeleLMeX_Diffusion.cpp | 109 ++++------------------------------ 1 file changed, 10 insertions(+), 99 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 0736819eb..721b4e85e 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -215,7 +215,6 @@ PeleLM::correctIsothermalBoundary( { BL_PROFILE("PeleLMeX::correctIsothermalBoundary()"); auto bcRecSpec = fetchBCRecArray(FIRSTSPEC, NUM_SPECIES); - bool need_explicit_fluxes = a_soretfluxes.empty(); Vector> soretfluxes(finest_level + 1); @@ -230,7 +229,7 @@ PeleLM::correctIsothermalBoundary( addSoretTerm( soretfluxes, soretfluxes, GetVecOfConstPtrs(getTempVect(a_time)), GetVecOfConstPtrs(getDiffusivityVect(a_time))); - } else { // have the lagged ones, alias to to them + } 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( @@ -295,98 +294,7 @@ PeleLM::correctIsothermalBoundary( } } } - - if (need_explicit_fluxes && m_use_wbar != 0) { // need to iterate to get the - // proper wbar flux - // check that we really don't have any wbar fluxes coming in - AMREX_ALWAYS_ASSERT(a_wbarfluxes.empty()); - Vector> wbarfluxes(finest_level + 1); - Vector residual(finest_level + 1); - const int iter_max = 10; - const Real tol = 1e-10; - for (int lev = 0; lev <= finest_level; lev++) { - // residual over all species - residual[lev].define(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); - for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { - wbarfluxes[lev][idim].define( - grids[lev], dmap[lev], NUM_SPECIES, 1, MFInfo(), Factory(lev)); - wbarfluxes[lev][idim].setVal(0.0); - } - } - for (int iter = 0; iter < iter_max; iter++) { - // based on existing gradient, compute the wbar flux - addWbarTerm( - GetVecOfArrOfPtrs(wbarfluxes), GetVecOfArrOfPtrs(wbarfluxes), - GetVecOfConstPtrs(getSpeciesVect(a_time)), - GetVecOfConstPtrs(getDensityVect(a_time)), - GetVecOfConstPtrs(getDiffusivityVect(a_time)), - GetVecOfConstPtrs(a_spec_boundary)); - for (int lev = 0; lev <= finest_level; ++lev) { - residual[lev].setVal(0.0); - auto* ldata_p = getLevelDataPtr(lev, a_time); - 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(*a_spec_boundary[lev], 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& rhoD_ec = beta_ec[idim].const_array(mfi); - auto const& flux_wbar = wbarfluxes[lev][idim].const_array(mfi); - auto const& flux_soret = soretfluxes[lev][idim]->const_array(mfi); - auto const& boundary_ar = a_spec_boundary[lev]->array(mfi); - auto const& residual_ar = residual[lev].array(mfi); - amrex::ParallelFor( - ebx, - [flux_soret, flux_wbar, rhoD_ec, boundary_ar, idim, edomain, - bc_lo, bc_hi, - residual_ar] 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++) { - // hold the old value - Real b_old = boundary_ar(idx[0], idx[1], idx[2], n); - // update the value - boundary_ar(idx[0], idx[1], idx[2], n) = - (flux_soret(i, j, k, n) + flux_wbar(i, j, k, n)) / - rhoD_ec(i, j, k, n); - // compute residual error - Real res = - std::abs(boundary_ar(idx[0], idx[1], idx[2], n) - b_old); - residual_ar(idx[0], idx[1], idx[2]) = - std::max(residual_ar(idx[0], idx[1], idx[2]), res); - } - } - }); - } - } - } - if (residual[finest_level].norm0(0, 1) < tol) { - break; - } - } - } + //TODO: wbar fluxes disabled for this case - boundary system becomes complex } void @@ -579,12 +487,14 @@ PeleLM::addWbarTerm( 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)); } + auto const* leosparm = eos_parms.device_parm(); const Box& domain = geom[lev].Domain(); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -603,14 +513,15 @@ PeleLM::addWbarTerm( gbx, [rho_arr, rhoY_arr, Wbar_arr, gradY_arr, Wbar_boundary_arr, domain, have_boundary, - phys_bc = m_phys_bc] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - getMwmixGivenRY(i, j, k, rho_arr, rhoY_arr, Wbar_arr); + 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++) { + 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 || @@ -619,9 +530,9 @@ PeleLM::addWbarTerm( 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); + getGradMwmixGivengradYMwmix(i, j, k, gradY_arr, Wbar_arr, Wbar_boundary_arr,leosparm); } } } From b52b30cd811443844d47e5b530582bf442cc649a Mon Sep 17 00:00:00 2001 From: Thomas Howarth Date: Fri, 15 Nov 2024 14:07:34 +0100 Subject: [PATCH 34/34] disable wbar --- Source/PeleLMeX_Diffusion.cpp | 21 +++++++++++---------- Source/PeleLMeX_K.H | 8 +++++--- Source/PeleLMeX_Setup.cpp | 8 +++++--- 3 files changed, 21 insertions(+), 16 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 93e29134e..5c43ced5e 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -414,7 +414,7 @@ PeleLM::correctIsothermalBoundary( } } } - //TODO: wbar fluxes disabled for this case - boundary system becomes complex + // TODO: wbar fluxes disabled for this case - boundary system becomes complex } void @@ -631,18 +631,18 @@ PeleLM::addWbarTerm( amrex::ParallelFor( - 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); + 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 + // 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++) { + 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 || @@ -651,9 +651,10 @@ PeleLM::addWbarTerm( 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); + getGradMwmixGivengradYMwmix( + i, j, k, gradY_arr, Wbar_arr, Wbar_boundary_arr, leosparm); } } } diff --git a/Source/PeleLMeX_K.H b/Source/PeleLMeX_K.H index b032bbb56..3ffbfa0ac 100644 --- a/Source/PeleLMeX_K.H +++ b/Source/PeleLMeX_K.H @@ -571,16 +571,18 @@ getGradMwmixGivengradYMwmix( int k, amrex::Array4 const& gradY, amrex::Array4 const& Mwmix, - amrex::Array4 const& gradMwmix) noexcept + amrex::Array4 const& gradMwmix, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + 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 + 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) diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 34ee0e6a0..d4c994154 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -103,9 +103,10 @@ PeleLM::Setup() << " 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" - << std::endl; + amrex::Print() + << " Imposing inhomogeneous Neumann conditions " + "for species on isothermal walls. WARNING: use_wbar disabled." + << std::endl; } } } else { @@ -437,6 +438,7 @@ PeleLM::readParameters() } 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