add advection/diffusion of a passive scalar

CMake/BuildROMSXExe.cmake

@@ -75,12 +75,11 @@ function(build_romsx_lib romsx_lib_name)
        ${SRC_DIR}/TimeIntegration/ROMSX_t3dmix.cpp
        ${SRC_DIR}/TimeIntegration/ROMSX_coriolis.cpp
        ${SRC_DIR}/TimeIntegration/ROMSX_prestep.cpp
-       ${SRC_DIR}/TimeIntegration/ROMSX_prestep_t_3d.cpp
-       ${SRC_DIR}/TimeIntegration/ROMSX_prestep_uv_3d.cpp
+       ${SRC_DIR}/TimeIntegration/ROMSX_prestep_t_advection.cpp
+       ${SRC_DIR}/TimeIntegration/ROMSX_prestep_diffusion.cpp
        ${SRC_DIR}/TimeIntegration/ROMSX_rhs_t_3d.cpp
        ${SRC_DIR}/TimeIntegration/ROMSX_rhs_uv_3d.cpp
        ${SRC_DIR}/TimeIntegration/ROMSX_rhs_uv_2d.cpp
-       ${SRC_DIR}/TimeIntegration/ROMSX_update_vel_3d.cpp
        ${SRC_DIR}/TimeIntegration/ROMSX_vert_visc_3d.cpp
        ${SRC_DIR}/TimeIntegration/ROMSX_set_massflux_3d.cpp
        ${SRC_DIR}/TimeIntegration/ROMSX_update_massflux_3d.cpp

Exec/Advection/CMakeLists.txt

@@ -0,0 +1,13 @@
+set(romsx_exe_name advection)
+
+add_executable(${romsx_exe_name} "")
+target_sources(${romsx_exe_name}
+   PRIVATE
+     prob.H
+     prob.cpp
+)
+
+target_include_directories(${romsx_exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
+
+include(${CMAKE_SOURCE_DIR}/CMake/BuildROMSXExe.cmake)
+build_romsx_exe(${romsx_exe_name})

Exec/Advection/GNUmakefile

@@ -0,0 +1,36 @@
+# AMReX
+COMP = gnu
+PRECISION = DOUBLE
+
+# 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
+
+TEST = TRUE
+USE_ASSERTION = TRUE
+
+# Debugging
+DEBUG = FALSE
+
+TEST = TRUE
+USE_ASSERTION = TRUE
+
+# GNU Make
+Bpack := ./Make.package
+Blocs := .
+ROMSX_HOME := ../..
+ROMSX_PROBLEM_DIR = $(ROMSX_HOME)/Exec/Advection
+USE_NETCDF = FALSE
+DEFINES += -DROMSX_USE_HISTORYFILE
+include $(ROMSX_HOME)/Exec/Make.ROMSX

Exec/Advection/Make.package

@@ -0,0 +1,2 @@
+CEXE_headers += prob.H
+CEXE_sources += prob.cpp

Exec/Advection/inputs

@@ -0,0 +1,67 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 10
+stop_time = 30000.0
+
+amrex.fpe_trap_invalid = 0
+
+fabarray.mfiter_tile_size = 16 16 1024
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_lo     =      0.     0.    -150.
+geometry.prob_hi     =  41000. 41000.       0.
+
+amr.n_cell           = 81     81      16
+
+# periodic in x to match WRF setup
+geometry.is_periodic = 1 1 0
+#ylo.type = "SlipWall"
+#yhi.type = "SlipWall"
+zlo.type = "SlipWall"
+zhi.type = "SlipWall"
+
+# TIME STEP CONTROL
+romsx.fixed_dt       = 300.0 # Timestep size (seconds)
+# NDTFAST  = 30.0 # Number of baratropic steps => 300.0/30.0 = 10.0
+romsx.fixed_fast_dt  = 10.0 # Baratropic timestep size (seconds)
+# romsx.fixed_fast_dt  = 300.0 # Baratropic timestep size (seconds) testing value
+romsx.fixed_ndtfast_ratio  = 30 # Baratropic timestep size (seconds)
+
+# DIAGNOSTICS & VERBOSITY
+romsx.sum_interval   = 10       # timesteps between computing mass
+romsx.v              = 0       # verbosity in ROMSX.cpp (0: none, 1: print boxes, etc, 2: print values)
+amr.v                = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+
+# CHECKPOINT FILES
+romsx.check_file      = chk        # root name of checkpoint file
+romsx.check_int       = -57600      # number of timesteps between checkpoints
+
+# PLOTFILES
+romsx.plot_file_1     = plt        # prefix of plotfile name
+romsx.plot_int_1      = 100          # number of timesteps between plotfiles
+romsx.plot_vars_1     = salt temp scalar x_velocity y_velocity z_velocity
+romsx.plotfile_type   = amrex
+
+# SOLVER CHOICE
+romsx.use_coriolis = false
+romsx.flat_bathymetry=true
+romsx.horizontal_advection_scheme = "upstream3" # upstream3 or centered4
+romsx.spatial_order = 2
+
+# Linear EOS parameters
+romsx.R0    = 1027.0  # background density value (Kg/m3) used in Linear Equation of State
+romsx.S0    = 35.0    # background salinity (nondimensional) constant
+romsx.T0    = 14.0    # background potential temperature (Celsius) constant
+romsx.Tcoef = 1.7e-4  # linear equation of state parameter (1/Celcius)
+romsx.Scoef = 0.0     # linear equation of state parameter (nondimensional)
+romsx.rho0  = 1025.0  # Mean density (Kg/m3) used when Boussinesq approx is inferred
+
+# Coriolis params
+romsx.coriolis_f0 = 0.0
+romsx.coriolis_beta = 0.0
+
+# PROBLEM PARAMETERS (velocity fields)
+prob.u_0   =1.0e-1
+prob.v_0   = -1.0e-1

Exec/Advection/prob.H

@@ -0,0 +1,16 @@
+#ifndef _PROB_H_
+#define _PROB_H_
+
+#include "AMReX_REAL.H"
+
+struct ProbParm {
+
+  amrex::Real u_0     = 0.0;
+  amrex::Real v_0     = 0.0;
+
+}; // namespace ProbParm
+
+extern ProbParm parms;
+
+#endif
+

Exec/Advection/prob.cpp

@@ -0,0 +1,208 @@
+#include "prob.H"
+#include "prob_common.H"
+
+#include "EOS.H"
+#include "AMReX_ParmParse.H"
+#include "AMReX_MultiFab.H"
+#include "IndexDefines.H"
+#include "DepthStretchTransform.H"
+
+using namespace amrex;
+
+ProbParm parms;
+
+void
+amrex_probinit(
+  const amrex_real* /*problo*/,
+  const amrex_real* /*probhi*/)
+{
+  // Parse params
+  ParmParse pp("prob");
+
+  pp.query("u_0", parms.u_0);
+  pp.query("v_0", parms.v_0);
+
+}
+
+/**
+ * \brief Initializes bathymetry h and surface height Zeta
+ */
+void
+init_custom_bathymetry (const Geometry& geom,
+                        MultiFab& mf_h,
+                        MultiFab& mf_Zt_avg1,
+                        const SolverChoice& m_solverChoice)
+{
+    const auto & geomdata = geom.data();
+    mf_h.setVal(geomdata.ProbHi(2));
+
+    const int Lm = geom.Domain().size()[0];
+    const int Mm = geom.Domain().size()[1];
+
+    //HACK HACK manually setting zeta to 0
+    mf_Zt_avg1.setVal(0.0);
+
+    for ( MFIter mfi(mf_h, TilingIfNotGPU()); mfi.isValid(); ++mfi )
+    {
+      Array4<Real> const& h  = (mf_h).array(mfi);
+
+      Box bx = mfi.tilebox();
+      Box gbx2 = bx;
+      gbx2.grow(IntVect(NGROW,NGROW,0));
+
+      bool NSPeriodic = geomdata.isPeriodic(1);
+      bool EWPeriodic = geomdata.isPeriodic(0);
+
+      Box gbx2D = gbx2;
+      gbx2D.makeSlab(2,0);
+
+      Gpu::streamSynchronize();
+      amrex::ParallelFor(gbx2,
+      [=] AMREX_GPU_DEVICE (int i, int j, int k)
+      {
+          h(i,j,0,0) = -geomdata.ProbLo(2);
+          if (k==0) {
+              h(i,j,0,1) = h(i,j,0,0);
+          }
+      });
+    } // mfi
+}
+
+void
+init_custom_prob(
+        const Box& bx,
+        Array4<Real      > const& state,
+        Array4<Real      > const& x_vel,
+        Array4<Real      > const& y_vel,
+        Array4<Real      > const& z_vel,
+        Array4<Real const> const& /*z_w*/,
+        Array4<Real const> const& z_r,
+        Array4<Real const> const& /*Hz*/,
+        Array4<Real const> const& /*h*/,
+        Array4<Real const> const& /*Zt_avg1*/,
+        GeometryData const& geomdata,
+        const SolverChoice& m_solverChoice)
+{
+    bool l_use_salt = m_solverChoice.use_salt;
+
+    const int khi = geomdata.Domain().bigEnd()[2];
+
+    AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);
+
+    ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    {
+        const auto prob_lo         = geomdata.ProbLo();
+        const auto prob_hi         = geomdata.ProbHi();
+        const auto dx              = geomdata.CellSize();
+
+        const Real z = z_r(i,j,k);
+
+        state(i, j, k, Temp_comp) = 1.;
+
+        state(i,j,k,Temp_comp)=m_solverChoice.T0; //+8.0*std::exp(z/50.0_rt);
+
+        // Set scalar = 0 everywhere
+        const Real xcent = 0.5*(prob_lo[0] + prob_hi[0]);
+        const Real ycent = 0.5*(prob_lo[1] + prob_hi[1]);
+
+        const Real x  = prob_lo[0] + (i + 0.5) * dx[0] - xcent;
+        const Real y  = prob_lo[1] + (j + 0.5) * dx[1] - ycent;
+        const Real r2 = x*x + y*y;
+        const Real rad = 0.1 * (prob_hi[0]-prob_lo[0]);
+        const Real radsq = rad*rad;
+
+        if (l_use_salt) {
+            state(i,j,k,Salt_comp)= m_solverChoice.S0;
+        }
+
+        state(i, j, k, Scalar_comp) = std::exp(-r2/(2.*radsq));
+    });
+
+  // Construct a box that is on x-faces
+  const Box& xbx = surroundingNodes(bx,0);
+  // Set the x-velocity
+  ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+  {
+        x_vel(i, j, k) = parms.u_0;
+  });
+
+  // Construct a box that is on y-faces
+  const Box& ybx = surroundingNodes(bx,1);
+
+  // Set the y-velocity
+  ParallelFor(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+  {
+        y_vel(i, j, k) = parms.v_0;
+  });
+
+  // Construct a box that is on z-faces
+  const Box& zbx = surroundingNodes(bx,2);
+
+  // Set the z-velocity
+  ParallelFor(zbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+  {
+      z_vel(i, j, k) = 0.0;
+  });
+
+  Gpu::streamSynchronize();
+}
+
+void
+init_custom_vmix(const Geometry& /*geom*/, MultiFab& mf_Akv, MultiFab& mf_Akt,
+                 MultiFab& mf_z_w, const SolverChoice& /*m_solverChoice*/)
+{
+    for ( MFIter mfi((mf_Akv), TilingIfNotGPU()); mfi.isValid(); ++mfi )
+    {
+      Array4<Real> const& Akv = (mf_Akv).array(mfi);
+      Array4<Real> const& Akt = (mf_Akt).array(mfi);
+      Array4<Real> const& z_w = (mf_z_w).array(mfi);
+      Box bx = mfi.tilebox();
+      bx.grow(IntVect(NGROW,NGROW,0));
+      Gpu::streamSynchronize();
+      amrex::ParallelFor(bx,
+      [=] AMREX_GPU_DEVICE (int i, int j, int k)
+      {
+        Akv(i,j,k) = 2.0e-03+8.0e-03*std::exp(z_w(i,j,k)/150.0);
+
+        Akt(i,j,k,Temp_comp) = 1.0e-6;
+        Akt(i,j,k,Salt_comp) = 1.0e-6;
+        Akt(i,j,k,Scalar_comp) = 0.0;
+      });
+    }
+}
+
+void
+init_custom_hmix(const Geometry& /*geom*/, MultiFab& mf_visc2_p, MultiFab& mf_visc2_r,
+                 MultiFab& mf_diff2, const SolverChoice& /*m_solverChoice*/)
+{
+    for ( MFIter mfi((mf_visc2_p), TilingIfNotGPU()); mfi.isValid(); ++mfi )
+    {
+      Array4<Real> const& visc2_p = (mf_visc2_p).array(mfi);
+      Array4<Real> const& visc2_r = (mf_visc2_r).array(mfi);
+      Array4<Real> const& diff2   = mf_diff2.array(mfi);
+      Box bx = mfi.tilebox();
+      bx.grow(IntVect(NGROW,NGROW,0));
+      Gpu::streamSynchronize();
+
+      int ncomp = mf_diff2.nComp();
+
+      amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+      {
+        visc2_p(i,j,k) = 5.0;
+        visc2_r(i,j,k) = 5.0;
+
+        for (int n = 0; n < ncomp; n++) {
+            diff2(i,j,k,n) = 0.0;
+        }
+      });
+    }
+}
+
+void
+init_custom_smflux(const Geometry& geom, const Real time, MultiFab& mf_sustr, MultiFab& mf_svstr,
+                   const SolverChoice& m_solverChoice)
+{
+    //It's possible these should be set to be nonzero only at the boundaries they affect
+    mf_sustr.setVal(0.0);
+    mf_svstr.setVal(0.0);
+}

Exec/CMakeLists.txt

@@ -4,3 +4,4 @@ include(${CMAKE_SOURCE_DIR}/CMake/BuildROMSXExe.cmake)
 build_romsx_lib(${romsx_lib_name})
 add_subdirectory(Seamount)
 add_subdirectory(Upwelling)
+add_subdirectory(Advection)

Exec/ParticlesOverSeaMount/prob.cpp

@@ -200,7 +200,10 @@ init_custom_vmix(const Geometry& /*geom*/, MultiFab& mf_Akv, MultiFab& mf_Akt, M
       ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
       {
           Akv(i,j,k) = 2.0e-03+8.0e-03*std::exp(z_w(i,j,k)/150.0);
-          Akt(i,j,k) = 1.0e-6_rt;
+
+          Akt(i,j,k,Temp_comp) = 1.0e-6_rt;
+          Akt(i,j,k,Salt_comp) = 1.0e-6_rt;
+          Akt(i,j,k,Scalar_comp) = 0.0_rt;
       });
     }
 }

Exec/Seamount/prob.cpp

@@ -165,7 +165,10 @@ init_custom_vmix(const Geometry& /*geom*/, MultiFab& mf_Akv, MultiFab& mf_Akt,
       amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
       {
         Akv(i,j,k) = 1.0e-5_rt;
-        Akt(i,j,k) = 1.0e-6_rt;
+
+        Akt(i,j,k,Temp_comp) = 1.0e-6_rt;
+        Akt(i,j,k,Salt_comp) = 1.0e-6_rt;
+        Akt(i,j,k,Scalar_comp) = 0.0_rt;
       });
     }
 }