Merge branch 'development' of github.com:erf-model/ERF into HEAD

CMake/BuildERFExe.cmake

@@ -36,8 +36,9 @@ function(build_erf_lib erf_lib_name)
 
   if(ERF_ENABLE_PARTICLES)
     target_sources(${erf_lib_name} PRIVATE
-                   ${SRC_DIR}/Particles/TracerPC.cpp
-                   ${SRC_DIR}/Particles/HydroPC.cpp)
+                   ${SRC_DIR}/Particles/ERFPCEvolve.cpp
+                   ${SRC_DIR}/Particles/ERFPCInitializations.cpp
+                   ${SRC_DIR}/Particles/ERFTracers.cpp)
     target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Particles)
     target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_PARTICLES)
   endif()

CMakeLists.txt

@@ -128,6 +128,11 @@ if(ERF_ENABLE_MPI)
   find_package(MPI REQUIRED)
 endif()
 
+configure_file(
+  ${CMAKE_SOURCE_DIR}/Source/ERF_Config.H.in
+  ${CMAKE_BINARY_DIR}/ERF_Config.H
+)
+
 # General information about machine, compiler, and build type
 message(STATUS "ERF Information:")
 message(STATUS "CMAKE_SYSTEM_NAME = ${CMAKE_SYSTEM_NAME}")

Source/ERF_Config.H.in

@@ -0,0 +1,15 @@
+
+#ifndef ERF_CONFIG_H_
+#define ERF_CONFIG_H_
+
+// define the radiation data directory
+#define RADIATION_DATA_DIR @ERF_RADIATION_DATA_DIR@
+
+// functions to get the radiation data directory
+inline
+std::string getRadiationDataDir()
+{
+   return "@ERF_RADIATION_DATA_DIR@";
+}
+
+#endif // ERF_CONFIG_H

Source/Particles/ERFPC.H

@@ -50,7 +50,6 @@ struct ERFParticlesAssignor
                                  amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
                                  const amrex::Box& domain ) const noexcept
     {
-        /* TODO: What is this about? */
         amrex::IntVect iv(
             AMREX_D_DECL(   int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])),
                             int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])),
@@ -61,6 +60,38 @@ struct ERFParticlesAssignor
     }
 };
 
+template <typename P>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void update_location_idata (P& p,
+                            amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
+                            amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
+                            const amrex::Array4<amrex::Real const>&  height_arr)
+{
+    amrex::IntVect iv( int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])),
+                       int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])),
+                       p.idata(ERFParticlesIntIdx::k) );
+
+    p.idata(ERFParticlesIntIdx::i) = iv[0];
+    p.idata(ERFParticlesIntIdx::j) = iv[1];
+
+    amrex::Real lx = (p.pos(0)-plo[0])*dxi[0] - static_cast<amrex::Real>(iv[0]);
+    amrex::Real ly = (p.pos(1)-plo[1])*dxi[1] - static_cast<amrex::Real>(iv[1]);
+    auto zlo = height_arr(iv[0]  ,iv[1]  ,iv[2]  ) * (1.0-lx) * (1.0-ly) +
+               height_arr(iv[0]+1,iv[1]  ,iv[2]  ) *      lx  * (1.0-ly) +
+               height_arr(iv[0]  ,iv[1]+1,iv[2]  ) * (1.0-lx) * ly +
+               height_arr(iv[0]+1,iv[1]+1,iv[2]  ) *      lx  * ly;
+    auto zhi = height_arr(iv[0]  ,iv[1]  ,iv[2]+1) * (1.0-lx) * (1.0-ly) +
+               height_arr(iv[0]+1,iv[1]  ,iv[2]+1) *      lx  * (1.0-ly) +
+               height_arr(iv[0]  ,iv[1]+1,iv[2]+1) * (1.0-lx) * ly +
+               height_arr(iv[0]+1,iv[1]+1,iv[2]+1) *      lx  * ly;
+
+    if (p.pos(2) > zhi) {
+        p.idata(ERFParticlesIntIdx::k) += 1;
+    } else if (p.pos(2) <= zlo) {
+        p.idata(ERFParticlesIntIdx::k) -= 1;
+    }
+}
+
 class ERFPC : public amrex::ParticleContainer<  ERFParticlesRealIdx::ncomps,
                                                 ERFParticlesIntIdx::ncomps,
                                                 0,
@@ -118,6 +149,17 @@ class ERFPC : public amrex::ParticleContainer<  ERFParticlesRealIdx::ncomps,
             return {AMREX_D_DECL("xvel","yvel","zvel"),"mass"};
         }
 
+        /*! Uses midpoint method to advance particles using flow velocity. */
+        virtual void AdvectWithFlow (  amrex::MultiFab*,
+                                       int,
+                                       amrex::Real,
+                                       const std::unique_ptr<amrex::MultiFab>& );
+
+        /*! Uses midpoint method to advance particles falling under gravity. */
+        virtual void AdvectWithGravity ( int,
+                                         amrex::Real,
+                                         const std::unique_ptr<amrex::MultiFab>& );
+
         /*! Specify if particles should advect with flow */
         inline void setAdvectWithFlow (bool a_flag)
         {
@@ -131,17 +173,6 @@ class ERFPC : public amrex::ParticleContainer<  ERFParticlesRealIdx::ncomps,
             m_advect_w_gravity = a_flag;
         }
 
-        /*! Uses midpoint method to advance particles using flow velocity. */
-        void AdvectWithFlow (    amrex::MultiFab*,
-                                 int,
-                                 amrex::Real,
-                                 const std::unique_ptr<amrex::MultiFab>& );
-
-        /*! Uses midpoint method to advance particles falling under gravity. */
-        void AdvectWithGravity ( int,
-                                 amrex::Real,
-                                 const std::unique_ptr<amrex::MultiFab>& );
-
     protected:
 
         bool m_advect_w_flow;               /*!< advect with flow velocity */

Source/Particles/ERFPCEvolve.cpp

@@ -103,14 +103,16 @@ void ERFPC::AdvectWithFlow ( MultiFab*                           a_umac,
                 ParticleType& p = p_pbox[i];
                 if (p.id() <= 0) { return; }
                 ParticleReal v[AMREX_SPACEDIM];
-                mac_interpolate(p, plo, dxi, umacarr, v);
+                mac_interpolate_mapped_z(p, plo, dxi, umacarr, zheight, v);
                 if (ipass == 0)
                 {
                     for (int dim=0; dim < AMREX_SPACEDIM; dim++)
                     {
                         p.rdata(dim) = p.pos(dim);
                         p.pos(dim) += static_cast<ParticleReal>(ParticleReal(0.5)*a_dt*v[dim]);
                     }
+                    // Update z-coordinate carried by the particle
+                    update_location_idata(p,plo,dxi,zheight);
                 }
                 else
                 {
@@ -119,35 +121,8 @@ void ERFPC::AdvectWithFlow ( MultiFab*                           a_umac,
                         p.pos(dim) = p.rdata(dim) + static_cast<ParticleReal>(a_dt*v[dim]);
                         p.rdata(dim) = v[dim];
                     }
-
                     // Update z-coordinate carried by the particle
-                    IntVect iv(
-                        AMREX_D_DECL(int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])),
-                                     int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])),
-                                     p.idata(0)));
-                    iv[0] += domain.smallEnd()[0];
-                    iv[1] += domain.smallEnd()[1];
-                    ParticleReal zlo, zhi;
-                    if (use_terrain) {
-                        Real lx = (p.pos(0)-plo[0])*dxi[0] - static_cast<Real>(iv[0]-domain.smallEnd()[0]);
-                        Real ly = (p.pos(1)-plo[1])*dxi[1] - static_cast<Real>(iv[1]-domain.smallEnd()[1]);
-                        zlo =  zheight(iv[0]  ,iv[1]  ,iv[2]  ) * (1.0-lx) * (1.0-ly) +
-                               zheight(iv[0]+1,iv[1]  ,iv[2]  ) *      lx  * (1.0-ly) +
-                               zheight(iv[0]  ,iv[1]+1,iv[2]  ) * (1.0-lx) * ly +
-                               zheight(iv[0]+1,iv[1]+1,iv[2]  ) *      lx  * ly;
-                        zhi =  zheight(iv[0]  ,iv[1]  ,iv[2]+1) * (1.0-lx) * (1.0-ly) +
-                               zheight(iv[0]+1,iv[1]  ,iv[2]+1) *      lx  * (1.0-ly) +
-                               zheight(iv[0]  ,iv[1]+1,iv[2]+1) * (1.0-lx) * ly +
-                               zheight(iv[0]+1,iv[1]+1,iv[2]+1) *      lx  * ly;
-                    } else {
-                        zlo =  iv[2]    * dx[2];
-                        zhi = (iv[2]+1) * dx[2];
-                    }
-                    if (p.pos(2) > zhi) { // need to be careful here
-                        p.idata(0) += 1;
-                    } else if (p.pos(2) <= zlo) {
-                        p.idata(0) -= 1;
-                    }
+                    update_location_idata(p,plo,dxi,zheight);
                 }
             });
         }
@@ -223,33 +198,7 @@ void ERFPC::AdvectWithGravity (  int                                 a_lev,
             p.rdata(ERFParticlesRealIdx::vz) -= (grav - drag) * half_dt;
 
             // also update z-coordinate here
-            IntVect iv(
-                AMREX_D_DECL(int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])),
-                             int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])),
-                             p.idata(0)));
-            iv[0] += domain.smallEnd()[0];
-            iv[1] += domain.smallEnd()[1];
-            ParticleReal zlo, zhi;
-            if (use_terrain) {
-                Real lx = (p.pos(0)-plo[0])*dxi[0] - static_cast<Real>(iv[0]-domain.smallEnd()[0]);
-                Real ly = (p.pos(1)-plo[1])*dxi[1] - static_cast<Real>(iv[1]-domain.smallEnd()[1]);
-                zlo =  zheight(iv[0]  ,iv[1]  ,iv[2]  ) * (1.0-lx) * (1.0-ly) +
-                       zheight(iv[0]+1,iv[1]  ,iv[2]  ) *      lx  * (1.0-ly) +
-                       zheight(iv[0]  ,iv[1]+1,iv[2]  ) * (1.0-lx) * ly +
-                       zheight(iv[0]+1,iv[1]+1,iv[2]  ) *      lx  * ly;
-                zhi =  zheight(iv[0]  ,iv[1]  ,iv[2]+1) * (1.0-lx) * (1.0-ly) +
-                       zheight(iv[0]+1,iv[1]  ,iv[2]+1) *      lx  * (1.0-ly) +
-                       zheight(iv[0]  ,iv[1]+1,iv[2]+1) * (1.0-lx) * ly +
-                       zheight(iv[0]+1,iv[1]+1,iv[2]+1) *      lx  * ly;
-            } else {
-                zlo =  iv[2]    * dx[2];
-                zhi = (iv[2]+1) * dx[2];
-            }
-            if (p.pos(2) > zhi) { // need to be careful here
-                p.idata(0) += 1;
-            } else if (p.pos(2) <= zlo) {
-                p.idata(0) -= 1;
-            }
+            update_location_idata(p,plo,dxi,zheight);
         });
     }