From e1b1bf4a66aa9fe6ae11ca6078749e70adf972e6 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 22 Aug 2023 19:00:56 -0700 Subject: [PATCH] rest of particle functionality -- but still need to replace z_phys_nd --- .../BoundaryConditions_bndryreg.cpp | 128 ------------------ Source/BoundaryConditions/ROMSX_FillPatch.cpp | 1 - Source/BoundaryConditions/ROMSX_PhysBCFunct.H | 5 - Source/IO/Checkpoint.cpp | 14 ++ Source/IO/Plotfile.cpp | 14 ++ Source/ROMSX.H | 6 +- Source/ROMSX.cpp | 24 +++- Source/TimeIntegration/ROMSX_Advance.cpp | 14 +- 8 files changed, 66 insertions(+), 140 deletions(-) delete mode 100644 Source/BoundaryConditions/BoundaryConditions_bndryreg.cpp diff --git a/Source/BoundaryConditions/BoundaryConditions_bndryreg.cpp b/Source/BoundaryConditions/BoundaryConditions_bndryreg.cpp deleted file mode 100644 index 1629b84c..00000000 --- a/Source/BoundaryConditions/BoundaryConditions_bndryreg.cpp +++ /dev/null @@ -1,128 +0,0 @@ -#include "ROMSX.H" - -using namespace amrex; - -// -// This routine uses data read in as BndryRegisters from a previous ERF run -// -// mf is the MultiFab to be filled -// time is the time at which the data should be filled -// bccomp is the index into both domain_bcs_type_bcr and bc_extdir_vals -// so this follows the BCVars enum -// -void -ROMSX::fill_from_bndryregs (const Vector& mfs, const Real time) -{ - // - // We now assume that if we read in on one face, we read in on all faces - // - AMREX_ALWAYS_ASSERT(m_r2d); - - int lev = 0; - const Box& domain = geom[lev].Domain(); - - const auto& dom_lo = amrex::lbound(domain); - const auto& dom_hi = amrex::ubound(domain); - - amrex::Vector>& bndry_data = m_r2d->interp_in_time(time); - - // xlo: ori = 0 - // ylo: ori = 1 - // zlo: ori = 2 - // xhi: ori = 3 - // yhi: ori = 4 - // zhi: ori = 5 - const auto& bdatxlo = (*bndry_data[0])[lev].const_array(); - const auto& bdatylo = (*bndry_data[1])[lev].const_array(); - const auto& bdatxhi = (*bndry_data[3])[lev].const_array(); - const auto& bdatyhi = (*bndry_data[4])[lev].const_array(); - - int bccomp; - - for (int var_idx = 0; var_idx < Vars::NumTypes; ++var_idx) - { - MultiFab& mf = *mfs[var_idx]; - const int icomp = 0; - const int ncomp = mf.nComp(); - - if (var_idx == Vars::xvel) { - bccomp = BCVars::xvel_bc; - } else if (var_idx == Vars::yvel) { - bccomp = BCVars::yvel_bc; - } else if (var_idx == Vars::zvel) { - bccomp = BCVars::zvel_bc; - } else if (var_idx == Vars::cons) { - bccomp = BCVars::cons_bc; - } - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(mf); mfi.isValid(); ++mfi) - { - const Array4& dest_arr = mf.array(mfi); - Box bx = mfi.validbox(); - - // x-faces - { - Box bx_xlo(bx); - bx_xlo.setSmall(0,dom_lo.x-1); bx_xlo.setBig(0,dom_lo.x-1); - bx_xlo.setSmall(1,dom_lo.y); bx_xlo.setBig(1,dom_hi.y); - bx_xlo.setSmall(2,dom_lo.z); bx_xlo.setBig(2,dom_hi.z); - if (var_idx == Vars::xvel) { - bx_xlo.setSmall(0,dom_lo.x); bx_xlo.setBig(0,dom_lo.x); - } else { - bx_xlo.setSmall(0,dom_lo.x-1); bx_xlo.setBig(0,dom_lo.x-1); - } - - Box bx_xhi(bx); - bx_xhi.setSmall(1,dom_lo.y ); bx_xhi.setBig(1,dom_hi.y ); - bx_xhi.setSmall(2,dom_lo.z ); bx_xhi.setBig(2,dom_hi.z ); - bx_xhi.setSmall(0,dom_hi.x+1); bx_xhi.setBig(0,dom_hi.x+1); - - ParallelFor( - bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - int jb = std::min(std::max(j,dom_lo.y),dom_hi.y); - int kb = std::min(std::max(k,dom_lo.z),dom_hi.z); - dest_arr(i,j,k,icomp+n) = bdatxlo(dom_lo.x-1,jb,kb,bccomp+n); - }, - bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - int jb = std::min(std::max(j,dom_lo.y),dom_hi.y); - int kb = std::min(std::max(k,dom_lo.z),dom_hi.z); - dest_arr(i,j,k,icomp+n) = bdatxhi(dom_hi.x+1,jb,kb,bccomp+n); - } - ); - } // x-faces - - // y-faces - { - Box bx_ylo(bx); - bx_ylo.setSmall(0,dom_lo.x); bx_ylo.setBig(0,dom_hi.x); - bx_ylo.setSmall(1,dom_lo.y); bx_ylo.setBig(1,dom_lo.y); - if (var_idx == Vars::yvel) { - bx_ylo.setSmall(1,dom_lo.y); bx_ylo.setBig(1,dom_lo.y); - } else { - bx_ylo.setSmall(1,dom_lo.y-1); bx_ylo.setBig(1,dom_lo.y-1); - } - - Box bx_yhi(bx); - bx_yhi.setSmall(0,dom_lo.x ); bx_yhi.setBig(0,dom_hi.x); - bx_yhi.setSmall(2,dom_lo.z ); bx_yhi.setBig(2,dom_hi.z); - bx_yhi.setSmall(1,dom_hi.y+1); bx_yhi.setBig(1,dom_hi.y+1); - - ParallelFor( - bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - int ib = std::min(std::max(i,dom_lo.x),dom_hi.x); - int kb = std::min(std::max(k,dom_lo.z),dom_hi.z); - dest_arr(i,j,k,icomp+n) = bdatylo(ib,dom_lo.y-1,kb,bccomp+n); - }, - bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - int ib = std::min(std::max(i,dom_lo.x),dom_hi.x); - int kb = std::min(std::max(k,dom_lo.z),dom_hi.z); - dest_arr(i,j,k,icomp+n) = bdatyhi(ib,dom_hi.y+1,kb,bccomp+n); - } - ); - } // y-faces - } // mf - } // var_idx -} diff --git a/Source/BoundaryConditions/ROMSX_FillPatch.cpp b/Source/BoundaryConditions/ROMSX_FillPatch.cpp index f5ef92fb..b18cf423 100644 --- a/Source/BoundaryConditions/ROMSX_FillPatch.cpp +++ b/Source/BoundaryConditions/ROMSX_FillPatch.cpp @@ -83,7 +83,6 @@ ROMSX::FillPatch (int lev, Real time, const Vector& mfs) (*physbcs[lev])(*mf,icomp_cons,ncomp_cons,ngvect_cons,time,cons_only); } /* - if (m_r2d) amrex::Abort("ReadBoundaryPlanes is not supported");//fill_from_bndryregs(mfs,time); #ifdef ROMSX_USE_NETCDF if (init_type == "real") amrex::Abort("This init type is not supported");//fill_from_wrfbdy(mfs,time); #endif diff --git a/Source/BoundaryConditions/ROMSX_PhysBCFunct.H b/Source/BoundaryConditions/ROMSX_PhysBCFunct.H index ceb34bb6..76aa57c9 100644 --- a/Source/BoundaryConditions/ROMSX_PhysBCFunct.H +++ b/Source/BoundaryConditions/ROMSX_PhysBCFunct.H @@ -99,11 +99,6 @@ public: const amrex::GpuArray dxInv, int icomp, int ncomp, amrex::Real time, int bccomp); - void fill_from_bndryregs (const int lev, const amrex::Box& bx, const amrex::Array4& dest_arr, - const int icomp, const int bccomp, const int ncomp, - const amrex::Box& domain, const amrex::BCRec* bc_ptr, - const amrex::Real time); - #ifdef ROMSX_USE_NETCDF void fill_from_wrfbdy (const int lev, const amrex::Box& bx, const amrex::Array4& dest_arr, const int icomp, const int bccomp, const int ncomp, diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp index 550b6364..403b499c 100644 --- a/Source/IO/Checkpoint.cpp +++ b/Source/IO/Checkpoint.cpp @@ -167,6 +167,12 @@ ROMSX::WriteCheckpointFile () const VisMF::Write(*(vec_zeta[lev]), amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "zeta")); VisMF::Write(*(vec_Zt_avg1[lev]), amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Zt_avg1")); } + +#ifdef ROMSX_USE_PARTICLES + if (use_tracer_particles) { + tracer_particles->Checkpoint(checkpointname, "tracers", true, tracer_particle_varnames); + } +#endif } void @@ -358,4 +364,12 @@ ROMSX::ReadCheckpointFile () vec_zeta[lev]->FillBoundary(geom[lev].periodicity()); vec_Zt_avg1[lev]->FillBoundary(geom[lev].periodicity()); } + +#ifdef ROMSX_USE_PARTICLES + if (use_tracer_particles) { + tracer_particles = std::make_unique(Geom(0), dmap[0], grids[0]); + std::string tracer_file("tracers"); + tracer_particles->Restart(restart_chkfile, tracer_file); + } +#endif } diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index e3e0d35b..9850dd9a 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -212,6 +212,13 @@ ROMSX::WritePlotFile (int which, Vector plot_var_names) varnames, Geom(), t_new[0], istep, refRatio()); writeJobInfo(plotfilename); + +#ifdef ROMSX_USE_PARTICLES + if (use_tracer_particles) { + tracer_particles->Checkpoint(plotfilename, "tracers", true, tracer_particle_varnames); + } +#endif + #ifdef ROMSX_USE_HDF5 } else if (plotfile_type == "hdf5" || plotfile_type == "HDF5") { amrex::Print() << "Writing plotfile " << plotfilename+"d01.h5" << "\n"; @@ -283,6 +290,13 @@ ROMSX::WritePlotFile (int which, Vector plot_var_names) WriteMultiLevelPlotfile(plotfilename, finest_level+1, GetVecOfConstPtrs(mf2), varnames, g2, t_new[0], istep, rr); writeJobInfo(plotfilename); + +#ifdef ROMSX_USE_PARTICLES + if (use_tracer_particles) { + tracer_particles->Checkpoint(plotfilename, "tracers", true, tracer_particle_varnames); + } +#endif + #ifdef ROMSX_USE_NETCDF } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") { for (int lev = 0; lev <= finest_level; ++lev) { diff --git a/Source/ROMSX.H b/Source/ROMSX.H index b0475c33..1a79035f 100644 --- a/Source/ROMSX.H +++ b/Source/ROMSX.H @@ -33,7 +33,7 @@ #include #include -#ifdef ERF_USE_PARTICLES +#ifdef ROMSX_USE_PARTICLES #include "TerrainFittedPC.H" #endif @@ -59,7 +59,7 @@ namespace InterpType { }; } -#ifdef ERF_USE_PARTICLES +#ifdef ROMSX_USE_PARTICLES typedef amrex::ParticleContainer TracerPC; #endif @@ -654,7 +654,7 @@ public: void set_weights (int lev); -#ifdef ERF_USE_PARTICLES +#ifdef ROMSX_USE_PARTICLES std::unique_ptr tracer_particles; static amrex::Vector tracer_particle_varnames; #endif diff --git a/Source/ROMSX.cpp b/Source/ROMSX.cpp index de6760aa..43822d4c 100644 --- a/Source/ROMSX.cpp +++ b/Source/ROMSX.cpp @@ -31,7 +31,6 @@ amrex::Real ROMSX::init_shrink = 1.0; amrex::Real ROMSX::change_max = 1.1; int ROMSX::fixed_ndtfast_ratio = 0; - // Type of mesh refinement algorithm int ROMSX::do_reflux = 0; int ROMSX::do_avg_down = 0; @@ -49,9 +48,16 @@ std::string ROMSX::plotfile_type = "amrex"; // init_type: "custom", "ideal", "real" std::string ROMSX::init_type = "custom"; +#ifdef ROMSX_USE_PARTICLES +bool ROMSX::use_tracer_particles = false; +amrex::Vector ROMSX::tracer_particle_varnames = {AMREX_D_DECL("xvel", "yvel", "zvel")}; +#endif + +#ifdef ROMSX_USE_NETCDF // NetCDF wrfinput (initialization) file(s) amrex::Vector> ROMSX::nc_init_file = {{""}}; // Must provide via input std::string ROMSX::nc_bdy_file; +#endif amrex::Vector BCNames = {"xlo", "ylo", "zlo", "xhi", "yhi", "zhi"}; @@ -274,6 +280,17 @@ ROMSX::InitData () AverageDown(); +#ifdef ROMSX_USE_PARTICLES + // Initialize tracer particles if required + if (use_tracer_particles) { + tracer_particles = std::make_unique(Geom(0), dmap[0], grids[0]); + + tracer_particles->InitParticles(*z_phys_nd[0]); + + Print() << "Initialized " << tracer_particles->TotalNumberOfParticles() << " tracer particles." << std::endl; + } +#endif + } else { // Restart from a checkpoint restart(); @@ -766,6 +783,11 @@ ROMSX::ReadParameters () pp.query("plot_file_2", plot_file_2); pp.query("plot_int_1", plot_int_1); pp.query("plot_int_2", plot_int_2); + +#ifdef ROMSX_USE_PARTICLES + // Tracer particle toggle + pp.query("use_tracer_particles", use_tracer_particles); +#endif } #ifdef ROMSX_USE_MULTIBLOCK diff --git a/Source/TimeIntegration/ROMSX_Advance.cpp b/Source/TimeIntegration/ROMSX_Advance.cpp index 0b5fddab..0152ed8e 100644 --- a/Source/TimeIntegration/ROMSX_Advance.cpp +++ b/Source/TimeIntegration/ROMSX_Advance.cpp @@ -612,23 +612,33 @@ ROMSX::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle } U_new.FillBoundary(geom[lev].periodicity()); - V_new.FillBoundary(geom[lev].periodicity()); U_old.FillBoundary(geom[lev].periodicity()); - V_old.FillBoundary(geom[lev].periodicity()); mf_temp.FillBoundary(geom[lev].periodicity()); mf_salt.FillBoundary(geom[lev].periodicity()); + mf_tempold.FillBoundary(geom[lev].periodicity()); mf_saltold.FillBoundary(geom[lev].periodicity()); vec_t3[lev]->FillBoundary(geom[lev].periodicity()); vec_s3[lev]->FillBoundary(geom[lev].periodicity()); + for (int lev = 0; lev <= finest_level; ++lev) { FillPatch(lev, t_new[lev], vars_new[lev]); } + +#ifdef ROMSX_USE_PARTICLES + // Update tracer particles on level 0 + if (lev == 0 && use_tracer_particles) { + MultiFab* Umac = {U_new, V_new, W_new}; + tracer_particles->AdvectWithUmac(Umac, lev, dt_lev, *z_phys_nd[0]); + } +#endif + + //We are not storing computed W aka Omega // MultiFab::Copy(W_new,mf_w,0,0,W_new.nComp(),IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0))); // W_new.FillBoundary();