Merge pull request #16 from seahorce-scidac/fix_for_particles

rest of particle functionality -- but still need to replace z_phys_nd

Source/BoundaryConditions/ROMSX_FillPatch.cpp

@@ -83,7 +83,6 @@ ROMSX::FillPatch (int lev, Real time, const Vector<MultiFab*>& 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

Source/BoundaryConditions/ROMSX_PhysBCFunct.H

@@ -99,11 +99,6 @@ public:
                           const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> dxInv,
                           int icomp, int ncomp, amrex::Real time, int bccomp);
 
-    void fill_from_bndryregs (const int lev,   const amrex::Box& bx, const amrex::Array4<amrex::Real>& 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<amrex::Real>& dest_arr,
                            const int icomp, const int bccomp, const int ncomp,

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<TerrainFittedPC>(Geom(0), dmap[0], grids[0]);
+       std::string tracer_file("tracers");
+       tracer_particles->Restart(restart_chkfile, tracer_file);
+   }
+#endif
 }

Source/IO/Plotfile.cpp

@@ -212,6 +212,13 @@ ROMSX::WritePlotFile (int which, Vector<std::string> 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<std::string> 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) {

Source/ROMSX.H

@@ -33,7 +33,7 @@
 #include <DataStruct.H>
 #include <Derive.H>
 
-#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<AMREX_SPACEDIM, 0, 0, 0> TracerPC;
 #endif
 
@@ -654,7 +654,7 @@ public:
 
     void set_weights (int lev);
 
-#ifdef ERF_USE_PARTICLES
+#ifdef ROMSX_USE_PARTICLES
     std::unique_ptr<TerrainFittedPC> tracer_particles;
     static amrex::Vector<std::string> tracer_particle_varnames;
 #endif

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<std::string> ROMSX::tracer_particle_varnames = {AMREX_D_DECL("xvel", "yvel", "zvel")};
+#endif
+
+#ifdef ROMSX_USE_NETCDF
 // NetCDF wrfinput (initialization) file(s)
 amrex::Vector<amrex::Vector<std::string>> ROMSX::nc_init_file = {{""}}; // Must provide via input
 std::string ROMSX::nc_bdy_file;
+#endif
 
 amrex::Vector<std::string> 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<TerrainFittedPC>(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

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();