From e63a23370e3f439781ea5602dca197ede5a39237 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Thu, 19 Dec 2024 15:47:57 -0700 Subject: [PATCH] RANS Part 1: Poisson wall distance (#2031) * Rename init_immersed_body to avoid confusion * Create wall dist array * Add RANS options, begin implementing walldist calc, add plotfile output Poisson solve and dist calc implemented for uniform dz for now * Fix boxarray, force 3D solve, call setLevelBC * Add masking to walldist poisson solve * Simplify code, use nodal solver instead * Update mask for nodal solver * Update walldist calc to properly handle nodal values Don't use poisson dist for thin bodies; directly calculate wall dist in trivial case of flat bottom boundary * Add clarifying print statements * Fix direct calc * GPU fix * Remove debug print, additional GPU fixes --------- Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> --- CMake/BuildERFExe.cmake | 1 + CMake/SetAmrexOptions.cmake | 2 +- Source/DataStructs/ERF_AdvStruct.H | 4 +- Source/DataStructs/ERF_TurbStruct.H | 9 + Source/ERF.H | 12 +- Source/ERF.cpp | 7 +- Source/ERF_MakeNewArrays.cpp | 30 +- Source/ERF_MakeNewLevel.cpp | 16 +- Source/IO/ERF_Plotfile.cpp | 4 + Source/LinearSolvers/ERF_PoissonWallDist.cpp | 310 ++++++++++++++++++ Source/LinearSolvers/Make.package | 1 + .../TimeIntegration/ERF_ComputeTimestep.cpp | 1 - 12 files changed, 375 insertions(+), 22 deletions(-) create mode 100644 Source/LinearSolvers/ERF_PoissonWallDist.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index dc4888ff7..113125f33 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -190,6 +190,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Utils/ERF_MomentumToVelocity.cpp ${SRC_DIR}/LinearSolvers/ERF_PoissonSolve.cpp ${SRC_DIR}/LinearSolvers/ERF_PoissonSolve_tb.cpp + ${SRC_DIR}/LinearSolvers/ERF_PoissonWallDist.cpp ${SRC_DIR}/LinearSolvers/ERF_ComputeDivergence.cpp ${SRC_DIR}/LinearSolvers/ERF_SolveWithGMRES.cpp ${SRC_DIR}/LinearSolvers/ERF_SolveWithMLMG.cpp diff --git a/CMake/SetAmrexOptions.cmake b/CMake/SetAmrexOptions.cmake index 038e3f0ab..ad1c504dc 100644 --- a/CMake/SetAmrexOptions.cmake +++ b/CMake/SetAmrexOptions.cmake @@ -29,7 +29,7 @@ set(AMReX_FORTRAN OFF) set(AMReX_LINEAR_SOLVERS ON) set(AMReX_LINEAR_SOLVERS_EM OFF) -set(AMReX_LINEAR_SOLVERS_INCFLO OFF) +set(AMReX_LINEAR_SOLVERS_INCFLO ON) set(AMReX_PARTICLES OFF) if(ERF_ENABLE_PARTICLES) diff --git a/Source/DataStructs/ERF_AdvStruct.H b/Source/DataStructs/ERF_AdvStruct.H index 5ce564338..9293cbc40 100644 --- a/Source/DataStructs/ERF_AdvStruct.H +++ b/Source/DataStructs/ERF_AdvStruct.H @@ -184,6 +184,7 @@ struct AdvChoice { pp.queryarr("zero_xflux_faces", zero_xflux); pp.queryarr("zero_yflux_faces", zero_yflux); pp.queryarr("zero_zflux_faces", zero_zflux); + have_zero_flux_faces = ((zero_xflux.size() > 0) || (zero_yflux.size() > 0) || (zero_zflux.size() > 0)); } void display() @@ -295,9 +296,10 @@ struct AdvChoice { amrex::Real moistscal_horiz_upw_frac = 1.0; amrex::Real moistscal_vert_upw_frac = 1.0; - // Thin immersed body + // Thin immersed bodies amrex::Vector zero_xflux; amrex::Vector zero_yflux; amrex::Vector zero_zflux; + bool have_zero_flux_faces = false; }; #endif diff --git a/Source/DataStructs/ERF_TurbStruct.H b/Source/DataStructs/ERF_TurbStruct.H index d649d61e7..dc9def089 100644 --- a/Source/DataStructs/ERF_TurbStruct.H +++ b/Source/DataStructs/ERF_TurbStruct.H @@ -5,6 +5,8 @@ AMREX_ENUM(LESType, None, Smagorinsky, Deardorff); +AMREX_ENUM(RANSType, None, kEqn); + AMREX_ENUM(PBLType, None, MYNN25, YSU); template @@ -36,6 +38,10 @@ struct TurbChoice { std::string les_type_string = "None"; query_one_or_per_level(pp, "les_type", les_type, lev, max_level); + // Which RANS closure? + std::string rans_type_string = "None"; + query_one_or_per_level(pp, "rans_type", rans_type, lev, max_level); + // Which PBL Closure static std::string pbl_type_string = "None"; query_one_or_per_level(pp, "pbl_type", pbl_type, lev, max_level); @@ -188,6 +194,9 @@ struct TurbChoice { amrex::Real theta_ref = 300.0; + // RANS type + RANSType rans_type; + // PBL model PBLType pbl_type; diff --git a/Source/ERF.H b/Source/ERF.H index 2458c9e54..6df3960ae 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -177,6 +177,9 @@ public: void project_velocities_tb (int lev, amrex::Real dt, amrex::Vector& vars, amrex::MultiFab& p); + // Calculate wall distance by solving a Poisson equation + void poisson_wall_dist (int lev); + #ifdef ERF_USE_FFT void solve_with_fft (int lev, amrex::MultiFab& rhs, amrex::MultiFab& p, amrex::Array& fluxes); @@ -404,7 +407,7 @@ public: void init_from_hse (int lev); - void init_immersed_body (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm); + void init_thin_body (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm); #ifdef ERF_USE_MULTIBLOCK // constructor used when ERF is created by a multiblock driver @@ -853,6 +856,9 @@ private: amrex::Vector> z_t_rk; + // Wall distance function + amrex::Vector> walldist; + // Map scale factors amrex::Vector> mapfac_m; amrex::Vector> mapfac_u; @@ -980,8 +986,10 @@ private: "Kmv","Kmh", // eddy diffusivity of heat "Khv","Khh", - // mynn pbl lengthscale + // turbulence lengthscale "Lturb", + // wall distance + "walldist", // moisture vars "qt", "qv", "qc", "qi", "qp", "qrain", "qsnow", "qgraup", "qsat", "rain_accum", "snow_accum", "graup_accum", diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 6b33761c0..f7a5acb4c 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -288,6 +288,9 @@ ERF::ERF_shared () z_t_rk.resize(nlevs_max); + // Wall distance + walldist.resize(nlevs_max); + // Mapfactors mapfac_m.resize(nlevs_max); mapfac_u.resize(nlevs_max); @@ -669,9 +672,7 @@ ERF::InitData_post () AverageDown(); } - if ((solverChoice.advChoice.zero_xflux.size() > 0) || - (solverChoice.advChoice.zero_yflux.size() > 0) || - (solverChoice.advChoice.zero_zflux.size() > 0)) + if (solverChoice.advChoice.have_zero_flux_faces) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 0, "Thin immersed body with refinement not currently supported."); diff --git a/Source/ERF_MakeNewArrays.cpp b/Source/ERF_MakeNewArrays.cpp index 791a18a80..5c3d91ad8 100644 --- a/Source/ERF_MakeNewArrays.cpp +++ b/Source/ERF_MakeNewArrays.cpp @@ -87,16 +87,26 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, z_t_rk[lev] = nullptr; } - // We use these area arrays regardless of terrain, EB or none of the above - detJ_cc[lev] = std::make_unique(ba,dm,1,1); - ax[lev] = std::make_unique(convert(ba,IntVect(1,0,0)),dm,1,1); - ay[lev] = std::make_unique(convert(ba,IntVect(0,1,0)),dm,1,1); - az[lev] = std::make_unique(convert(ba,IntVect(0,0,1)),dm,1,1); + // We use these area arrays regardless of terrain, EB or none of the above + detJ_cc[lev] = std::make_unique(ba,dm,1,1); + ax[lev] = std::make_unique(convert(ba,IntVect(1,0,0)),dm,1,1); + ay[lev] = std::make_unique(convert(ba,IntVect(0,1,0)),dm,1,1); + az[lev] = std::make_unique(convert(ba,IntVect(0,0,1)),dm,1,1); - detJ_cc[lev]->setVal(1.0); - ax[lev]->setVal(1.0); - ay[lev]->setVal(1.0); - az[lev]->setVal(1.0); + detJ_cc[lev]->setVal(1.0); + ax[lev]->setVal(1.0); + ay[lev]->setVal(1.0); + az[lev]->setVal(1.0); + + // ******************************************************************************************** + // Create wall distance array for RANS modeling + // ******************************************************************************************** + if (solverChoice.turbChoice[lev].rans_type != RANSType::None) { + walldist[lev] = std::make_unique(ba,dm,1,1); + walldist[lev]->setVal(1e23); + } else { + walldist[lev] = nullptr; + } // ******************************************************************************************** // These are the persistent containers for the old and new data @@ -237,7 +247,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, #if defined(ERF_USE_WINDFARM) //********************************************************* - // Variables for Ftich model for windfarm parametrization + // Variables for Fitch model for windfarm parametrization //********************************************************* if (solverChoice.windfarm_type == WindFarmType::Fitch){ vars_windfarm[lev].define(ba, dm, 5, ngrow_state); // V, dVabsdt, dudt, dvdt, dTKEdt diff --git a/Source/ERF_MakeNewLevel.cpp b/Source/ERF_MakeNewLevel.cpp index fdfdf6eb5..addc31964 100644 --- a/Source/ERF_MakeNewLevel.cpp +++ b/Source/ERF_MakeNewLevel.cpp @@ -101,7 +101,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, // ******************************************************************************************** // Thin immersed body // ******************************************************************************************* - init_immersed_body(lev, ba, dm); + init_thin_body(lev, ba, dm); // ******************************************************************************************** // Initialize the integrator class @@ -144,6 +144,14 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, if (solverChoice.do_forest_drag) { m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } if (solverChoice.do_terrain_drag) { m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + + //******************************************************************************************** + // Create wall distance field for RANS model (depends upon z_phys) + // ******************************************************************************************* + if (solverChoice.turbChoice[lev].rans_type != RANSType::None) { + poisson_wall_dist(lev); + } + //******************************************************************************************** // Microphysics // ******************************************************************************************* @@ -496,7 +504,7 @@ ERF::ClearLevel (int lev) } void -ERF::init_immersed_body (int lev, const BoxArray& ba, const DistributionMapping& dm) +ERF::init_thin_body (int lev, const BoxArray& ba, const DistributionMapping& dm) { //******************************************************************************************** // Thin immersed body @@ -544,7 +552,7 @@ ERF::init_immersed_body (int lev, const BoxArray& ba, const DistributionMapping& } if (solverChoice.advChoice.zero_yflux.size() > 0) { - amrex::Print() << "Setting up thin interface boundary for " + amrex::Print() << "Setting up thin immersed body for " << solverChoice.advChoice.zero_yflux.size() << " yfaces" << std::endl; BoxArray ba_yf(ba); ba_yf.surroundingNodes(1); @@ -576,7 +584,7 @@ ERF::init_immersed_body (int lev, const BoxArray& ba, const DistributionMapping& } if (solverChoice.advChoice.zero_zflux.size() > 0) { - amrex::Print() << "Setting up thin interface boundary for " + amrex::Print() << "Setting up thin immersed body for " << solverChoice.advChoice.zero_zflux.size() << " zfaces" << std::endl; BoxArray ba_zf(ba); ba_zf.surroundingNodes(2); diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index c54c61a1f..5ce963f46 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -1045,6 +1045,10 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Turb_lengthscale,mf_comp,1,0); mf_comp ++; } + if (containerHasElement(plot_var_names, "walldist")) { + MultiFab::Copy(mf[lev],*walldist[lev],0,mf_comp,1,0); + mf_comp ++; + } // TODO: The size of the q variables can vary with different // moisture models. Therefore, certain components may diff --git a/Source/LinearSolvers/ERF_PoissonWallDist.cpp b/Source/LinearSolvers/ERF_PoissonWallDist.cpp new file mode 100644 index 000000000..2555e9204 --- /dev/null +++ b/Source/LinearSolvers/ERF_PoissonWallDist.cpp @@ -0,0 +1,310 @@ +#include "ERF.H" +#include "ERF_Utils.H" + +#include +#include + +using namespace amrex; + +/** + * Calculate wall distances using the Poisson equation + * + * See Tucker, P. G. (2003). Differential equation-based wall distance + * computation for DES and RANS. Journal of Computational Physics, + * 190(1), 229–248. https://doi.org/10.1016/S0021-9991(03)00272-9 + */ +void ERF::poisson_wall_dist (int lev) +{ + BL_PROFILE("ERF::poisson_wall_dist()"); + + auto const& geomdata = geom[lev]; + + if (solverChoice.mesh_type == MeshType::ConstantDz) { +// Comment this out to test the wall dist calc in the trivial case: +//#if 0 + Print() << "Directly calculating direct wall distance for constant dz" << std::endl; + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + for (MFIter mfi(*walldist[lev]); mfi.isValid(); ++mfi) { + const Box& bx = mfi.validbox(); + auto dist_arr = walldist[lev]->array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + dist_arr(i, j, k) = prob_lo[2] + (k + 0.5) * dx[2]; + }); + } + return; +//#endif + } else if (solverChoice.mesh_type == MeshType::StretchedDz) { + // TODO: Handle this trivial case + Error("Wall dist calc not implemented with grid stretching yet"); + } else { + // TODO + Error("Wall dist calc not implemented over terrain yet"); + } + + Print() << "Calculating Poisson wall distance" << std::endl; + + // Make sure the solver only sees the levels over which we are solving + BoxArray nba = walldist[lev]->boxArray(); + nba.surroundingNodes(); + Vector geom_tmp; geom_tmp.push_back(geom[lev]); + Vector ba_tmp; ba_tmp.push_back(nba); + Vector dm_tmp; dm_tmp.push_back(walldist[lev]->DistributionMap()); + + Vector rhs; + Vector phi; + +#ifdef ERF_USE_EB + Error("Wall dist calc not implemented for EB"; +#else + rhs.resize(1); rhs[0].define(ba_tmp[0], dm_tmp[0], 1, 0); + phi.resize(1); phi[0].define(ba_tmp[0], dm_tmp[0], 1, 1); +#endif + + rhs[0].setVal(-1.0); + + // Define an overset mask to set dirichlet nodes on walls + iMultiFab mask(ba_tmp[0], dm_tmp[0], 1, 0); + Vector overset_mask = {&mask}; + + auto const dom_lo = lbound(geom[lev].Domain()); + auto const dom_hi = ubound(geom[lev].Domain()); + + // **************************************************************************** + // Initialize phi + // (It is essential that we do this in order to fill the corners; this is + // used if we include blanking.) + // **************************************************************************** + phi[0].setVal(0.0); + + // **************************************************************************** + // Interior boundaries are marked with phi=0 + // **************************************************************************** + // Overset mask is 0/1: 1 means the node is an unknown. 0 means it's known. + mask.setVal(1); + if (solverChoice.advChoice.have_zero_flux_faces) { + Warning("Poisson distance is inaccurate for bodies in open domains that are small compared to the domain size, skipping..."); + walldist[lev]->setVal(1e34); + return; +#if 0 + Gpu::DeviceVector xfacelist, yfacelist, zfacelist; + + xfacelist.resize(solverChoice.advChoice.zero_xflux.size()); + yfacelist.resize(solverChoice.advChoice.zero_yflux.size()); + zfacelist.resize(solverChoice.advChoice.zero_zflux.size()); + + if (xfacelist.size() > 0) { + Gpu::copy(amrex::Gpu::hostToDevice, + solverChoice.advChoice.zero_xflux.begin(), + solverChoice.advChoice.zero_xflux.end(), + xfacelist.begin()); + Print() << " masking interior xfaces" << std::endl; + } + if (yfacelist.size() > 0) { + Gpu::copy(amrex::Gpu::hostToDevice, + solverChoice.advChoice.zero_yflux.begin(), + solverChoice.advChoice.zero_yflux.end(), + yfacelist.begin()); + Print() << " masking interior yfaces" << std::endl; + } + if (zfacelist.size() > 0) { + Gpu::copy(amrex::Gpu::hostToDevice, + solverChoice.advChoice.zero_zflux.begin(), + solverChoice.advChoice.zero_zflux.end(), + zfacelist.begin()); + Print() << " masking interior zfaces" << std::endl; + } + + for (MFIter mfi(phi[0]); mfi.isValid(); ++mfi) { + const Box& bx = mfi.validbox(); + + auto phi_arr = phi[0].array(mfi); + auto mask_arr = mask.array(mfi); + + if (xfacelist.size() > 0) { + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + for (int iface=0; iface < xfacelist.size(); ++iface) { + if ((i == xfacelist[iface][0]) && + (j == xfacelist[iface][1]) && + (k == xfacelist[iface][2])) + { + mask_arr(i, j , k ) = 0; + mask_arr(i, j , k+1) = 0; + mask_arr(i, j+1, k ) = 0; + mask_arr(i, j+1, k+1) = 0; + } + } + }); + } + + if (yfacelist.size() > 0) { + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + for (int iface=0; iface < yfacelist.size(); ++iface) { + if ((i == yfacelist[iface][0]) && + (j == yfacelist[iface][1]) && + (k == yfacelist[iface][2])) + { + mask_arr(i , j, k ) = 0; + mask_arr(i , j, k+1) = 0; + mask_arr(i+1, j, k ) = 0; + mask_arr(i+1, j, k+1) = 0; + } + } + }); + } + + if (zfacelist.size() > 0) { + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + for (int iface=0; iface < zfacelist.size(); ++iface) { + if ((i == xfacelist[iface][0]) && + (j == xfacelist[iface][1]) && + (k == xfacelist[iface][2])) + { + mask_arr(i , j , k) = 0; + mask_arr(i , j+1, k) = 0; + mask_arr(i+1, j , k) = 0; + mask_arr(i+1, j+1, k) = 0; + } + } + }); + } + } +#endif + } + + // **************************************************************************** + // Setup BCs, with solid domain boundaries being dirichlet + // We assume that the zlo boundary corresponds to the land surface + // **************************************************************************** + amrex::Array bc3d_lo, bc3d_hi; + Orientation zlo(Direction::z, Orientation::low); + bool havewall{false}; + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + if (geom[0].isPeriodic(dir)) { + bc3d_lo[dir] = LinOpBCType::Periodic; + bc3d_hi[dir] = LinOpBCType::Periodic; + } else { + bc3d_lo[dir] = LinOpBCType::Neumann; + bc3d_hi[dir] = LinOpBCType::Neumann; + } + } + if ( ( phys_bc_type[zlo] == ERF_BC::MOST ) || + ( phys_bc_type[zlo] == ERF_BC::no_slip_wall ) )/*|| + ((phys_bc_type[zlo] == ERF_BC::slip_wall) && (dom_hi.z > dom_lo.z)) )*/ + { + Print() << " Poisson zlo BC is dirichlet" << std::endl; + bc3d_lo[2] = LinOpBCType::Dirichlet; + havewall = true; + } + Print() << " bc lo : " << bc3d_lo << std::endl; + Print() << " bc hi : " << bc3d_hi << std::endl; + + if (!solverChoice.advChoice.have_zero_flux_faces && !havewall) { + Error("No solid boundaries in the computational domain"); + } + + LPInfo info; +/* Nodal solver cannot have hidden dimensions */ +#if 0 + // Allow a hidden direction if the domain is one cell wide + if (dom_lo.x == dom_hi.x) { + info.setHiddenDirection(0); + Print() << " domain is 2D in yz" << std::endl; + } else if (dom_lo.y == dom_hi.y) { + info.setHiddenDirection(1); + Print() << " domain is 2D in xz" << std::endl; + } else if (dom_lo.z == dom_hi.z) { + info.setHiddenDirection(2); + Print() << " domain is 2D in xy" << std::endl; + } +#endif + + // **************************************************************************** + // Solve nodal masked Poisson problem with MLMG + // TODO: different solver for terrain? + // **************************************************************************** + const Real reltol = solverChoice.poisson_reltol; + const Real abstol = solverChoice.poisson_abstol; + + Real sigma = 1.0; + MLNodeLaplacian mlpoisson(geom_tmp, ba_tmp, dm_tmp, info, {}, sigma); + + mlpoisson.setDomainBC(bc3d_lo, bc3d_hi); + + if (lev > 0) { + mlpoisson.setCoarseFineBC(nullptr, ref_ratio[lev-1], LinOpBCType::Neumann); + } + + mlpoisson.setLevelBC(0, nullptr); + + mlpoisson.setOversetMask(0, mask); + + // Solve + MLMG mlmg(mlpoisson); + int max_iter = 100; + mlmg.setMaxIter(max_iter); + + mlmg.setVerbose(mg_verbose); + mlmg.setBottomVerbose(0); + + mlmg.solve(GetVecOfPtrs(phi), + GetVecOfConstPtrs(rhs), + reltol, abstol); + + // Now overwrite with periodic fill outside domain and fine-fine fill inside + phi[0].FillBoundary(geom[lev].periodicity()); + + // **************************************************************************** + // Compute grad(phi) to get distances + // - Note that phi is nodal and walldist is cell-centered + // - TODO: include terrain metrics for dphi/dz + // **************************************************************************** + for (MFIter mfi(*walldist[lev]); mfi.isValid(); ++mfi) { + const Box& bx = mfi.validbox(); + + const auto invCellSize = geomdata.InvCellSizeArray(); + + auto const& phi_arr = phi[0].const_array(mfi); + auto dist_arr = walldist[lev]->array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + Real dpdx{0}, dpdy{0}, dpdz{0}; + + // dphi/dx + if (dom_lo.x != dom_hi.x) { + dpdx = 0.25 * invCellSize[0] * ( + (phi_arr(i+1, j , k ) - phi_arr(i, j , k )) + + (phi_arr(i+1, j , k+1) - phi_arr(i, j , k+1)) + + (phi_arr(i+1, j+1, k ) - phi_arr(i, j+1, k )) + + (phi_arr(i+1, j+1, k+1) - phi_arr(i, j+1, k+1)) ); + } + + // dphi/dy + if (dom_lo.y != dom_hi.y) { + dpdy = 0.25 * invCellSize[1] * ( + (phi_arr(i , j+1, k ) - phi_arr(i , j, k )) + + (phi_arr(i , j+1, k+1) - phi_arr(i , j, k+1)) + + (phi_arr(i+1, j+1, k ) - phi_arr(i+1, j, k )) + + (phi_arr(i+1, j+1, k+1) - phi_arr(i+1, j, k+1)) ); + } + + // dphi/dz + if (dom_lo.z != dom_hi.z) { + dpdz = 0.25 * invCellSize[2] * ( + (phi_arr(i , j , k+1) - phi_arr(i , j , k)) + + (phi_arr(i , j+1, k+1) - phi_arr(i , j+1, k)) + + (phi_arr(i+1, j , k+1) - phi_arr(i+1, j , k)) + + (phi_arr(i+1, j+1, k+1) - phi_arr(i+1, j+1, k)) ); + } + + Real dp_dot_dp = dpdx*dpdx + dpdy*dpdy + dpdz*dpdz; + Real phi_avg = 0.125 * ( + phi_arr(i , j , k ) + phi_arr(i , j , k+1) + phi_arr(i , j+1, k ) + phi_arr(i , j+1, k+1) + + phi_arr(i+1, j , k ) + phi_arr(i+1, j , k+1) + phi_arr(i+1, j+1, k ) + phi_arr(i+1, j+1, k+1) ); + dist_arr(i, j, k) = -std::sqrt(dp_dot_dp) + std::sqrt(dp_dot_dp + 2*phi_avg); + + // DEBUG: output phi instead + //dist_arr(i, j, k) = phi_arr(i, j, k); + }); + } +} diff --git a/Source/LinearSolvers/Make.package b/Source/LinearSolvers/Make.package index e966146b8..a499733e3 100644 --- a/Source/LinearSolvers/Make.package +++ b/Source/LinearSolvers/Make.package @@ -5,6 +5,7 @@ CEXE_sources += ERF_TerrainPoisson.cpp CEXE_sources += ERF_ComputeDivergence.cpp CEXE_sources += ERF_PoissonSolve.cpp CEXE_sources += ERF_PoissonSolve_tb.cpp +CEXE_sources += ERF_PoissonWallDist.cpp CEXE_sources += ERF_SolveWithGMRES.cpp CEXE_sources += ERF_SolveWithMLMG.cpp diff --git a/Source/TimeIntegration/ERF_ComputeTimestep.cpp b/Source/TimeIntegration/ERF_ComputeTimestep.cpp index 5c79fb85a..bd43780be 100644 --- a/Source/TimeIntegration/ERF_ComputeTimestep.cpp +++ b/Source/TimeIntegration/ERF_ComputeTimestep.cpp @@ -204,7 +204,6 @@ ERF::estTimeStep (int level, long& dt_fast_ratio) const } else if (fixed_dt[level] > 0.) { // Max CFL_c = 1.0 for substeps by default, but we enforce a min of 4 substeps auto dt_sub_max = (estdt_comp/cfl * sub_cfl); - Print() << "fixed_dt="<