Moving bilinear interpolation to Utils

Exec/DevTests/Hurricane/ERF_Prob.cpp

@@ -1,6 +1,8 @@
 #include "ERF_Prob.H"
 #include <ERF_Constants.H>
 
+#include "ERF_Interpolation_Bilinear.H"
+
 using namespace amrex;
 
 std::unique_ptr<ProblemBase>
@@ -84,80 +86,6 @@ Problem::erf_init_dens_hse_moist (MultiFab& rho_hse,
     } // no terrain
 }
 
-AMREX_FORCE_INLINE
-AMREX_GPU_HOST_DEVICE
-int get_single_index(int i, int j, int k,
-                     int nx, int ny)
-{
-
-    int si = k*nx*ny + j*nx + i;
-    return si;
-}
-
-AMREX_FORCE_INLINE
-AMREX_GPU_HOST_DEVICE
-void bilinear_interpolation(const Real* xvec, const Real* yvec, const Real* zvec,
-                            const Real dxvec, const Real dyvec,
-                            const int nx, const int ny, const int nz,
-                            const Real x, const Real y, const Real z,
-                            const Real* varvec,
-                            Real& tmp_var)
-{
-    int iloc=-1, jloc=-1, kloc=-1;
-    for(int k=0;k<nz;k++){
-        if(zvec[k] > z){
-            kloc = k-1;
-            break;
-        }
-        else if (zvec[k] == z) {
-            kloc = k;
-        }
-
-    }
-        iloc = std::floor((x-xvec[0])/dxvec);
-        jloc = std::floor((y-yvec[0])/dyvec);
-
-        if(iloc > nx-1 or iloc < 0 or
-           jloc > ny-1 or iloc < 0 or
-           kloc > nz-1 or kloc < 0){
-            //std::cout << "The value of iloc, jloc, kloc is " << iloc << " " << jloc <<
-            //                                                    kloc << "\n";
-            //exit(0);
-        }
-
-        Real xlo = xvec[0] + iloc*dxvec;
-        Real ylo = yvec[0] + jloc*dyvec;
-        Real zlo = zvec[kloc];
-
-        /*std::cout << "dxvec and dyvec are " << dxvec << " " << dyvec << "\n";
-        std::cout << "The value of xvec0, yvec0 is " << xvec[0] << " " << yvec[0] <<  "\n";
-        std::cout << "The value of x-xvec0/dxvec, y-yvec0/dyvec is " << x-xvec[0] << " " << y-yvec[0] <<  "\n";
-        std::cout << "The value of x-xvec0/dxvec, y-yvec0/dyvec is " << (x-xvec[0])/dxvec << " " << (y-yvec[0])/dyvec <<  "\n";
-        std::cout << "The value of xlo, ylo, zlo is " << xlo << " " << ylo << " " <<  zlo << "\n";
-        std::cout << "The value of x, y, z is " << x << " " << y << " " <<  z << "\n";
-        std::cout << "iloc, jloc, kloc = " << iloc << " " << jloc << " " << kloc << "\n";*/
-
-        Real w_x = (x - xlo)/dxvec;
-        Real w_y = (y - ylo)/dyvec;
-        Real w_z = (z - zlo)/(zvec[kloc+1] - zvec[kloc]);
-
-        int ind0 = get_single_index(iloc,jloc,kloc,nx,ny);
-        int ind1 = get_single_index(iloc+1,jloc,kloc,nx,ny);
-        int ind2 = get_single_index(iloc+1,jloc+1,kloc,nx,ny);
-        int ind3 = get_single_index(iloc,jloc+1,kloc,nx,ny);
-        int ind4 = get_single_index(iloc,jloc,kloc+1,nx,ny);
-        int ind5 = get_single_index(iloc+1,jloc,kloc+1,nx,ny);
-        int ind6 = get_single_index(iloc+1,jloc+1,kloc+1,nx,ny);
-        int ind7 = get_single_index(iloc,jloc+1,kloc+1,nx,ny);
-
-        tmp_var = (1-w_x)*(1-w_y)*(1-w_z)*varvec[ind0] + w_x*(1-w_y)*(1-w_z)*varvec[ind1] +
-                  w_x*w_y*(1-w_z)*varvec[ind2] + (1-w_x)*w_y*(1-w_z)*varvec[ind3] +
-                  (1-w_x)*(1-w_y)*w_z*varvec[ind4] + w_x*(1-w_y)*w_z*varvec[ind5] +
-                  w_x*w_y*w_z*varvec[ind6] + (1-w_x)*w_y*w_z*varvec[ind7];
-
-        //std::cout << "Variable value is " << tmp_var << "\n";
-}
-
 void
 Problem::init_custom_pert (
     const Box& bx,

Exec/DevTests/Hurricane/README.md

@@ -8,20 +8,20 @@ This folder contains examples for hurricane simulations from real weather data.
 1. Download the initial condition file 
 
 ```
-wget https://zenodo.org/record/14880478/files/ERF_IC_Henri_gdas1.fnl0p25.2021081906.f00.bin
+wget https://zenodo.org/record/15043093/files/ERF_IC_HenriERA5_20210819_VeryLarge.bin
 ```
 
 2. `make -j8`
-3. Run with `inputs_20210819_Henri`
+3. Run with `inputs_20210819_Henri_NoAMR`
 
 ## Hurricane Laura
 
 1. Download the initial condition file 
 
 ```
-wget https://zenodo.org/record/14880478/files/ERF_IC_Laura_2020082600.bin
+wget https://zenodo.org/record/15043093/files/ERF_IC_Laura_LargeDomain.bin
 ```
 
 2. `make -j8`.
-3. Run with `inputs_20200826_Laura`
+3. Run with `inputs_20200826_Laura_NoAMR` or `inputs_20200826_Laura_2LevelAMR`.
 

Exec/DevTests/Hurricane/inputs_20210821_Henri → Exec/DevTests/Hurricane/inputs_20200826_Laura_2LevelAMR

@@ -8,11 +8,16 @@ amrex.fpe_trap_invalid = 1
 fabarray.mfiter_tile_size = 2048 1024 2048
 
 # PROBLEM SIZE & GEOMETRY
-geometry.prob_lo     = 1000001 2000001    0.
-geometry.prob_hi     = 4249999 4999999 25000.
-amr.n_cell           =  500    500    100    # dx=dy=dz=250 m
+geometry.prob_lo     = -1999999  1499999    0.
+geometry.prob_hi     =  3999999  5259999 25000.
+amr.n_cell           =  1000    500    100    # dx=dy=dz=250 m
 
-erf.IC_file = "ERF_IC_Henri_gdas1.fnl0p25.2021082100.f00.bin"
+amr.max_level = 2
+erf.regrid_int = 100
+amr.ref_ratio_vect = 2 2 1 2 2 1
+erf.coupling_type  = "TwoWay"
+
+erf.IC_file = "ERF_IC_Laura_LargeDomain.bin"
 
 # periodic in x to match WRF setup
 # - as an alternative, could use symmetry at x=0 and outflow at x=25600
@@ -24,6 +29,11 @@ yhi.type = "Open"
 zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
+erf.use_coriolis = true
+erf.coriolis_3d = true
+erf.latitude = 30.0
+erf.rotational_time_period = 86165.734375
+
 #erf.rayleigh_damp_T = true
 #erf.rayleigh_zdamp = 5000.0
 #erf.rayleigh_dampcoef = 0.005
@@ -34,12 +44,9 @@ erf.fixed_fast_dt  = 2.5     # fixed time step [s]
 
 # DIAGNOSTICS & VERBOSITY
 erf.sum_interval   = 1       # timesteps between computing mass
-erf.v              = 1       # verbosity in ERF.cpp
+erf.v              = 0       # verbosity in ERF.cpp
 amr.v              = 1       # verbosity in Amr.cpp
 
-# REFINEMENT / REGRIDDING
-amr.max_level       = 0       # maximum level number allowed
-
 # CHECKPOINT FILES
 amr.check_file      = chk        # root name of checkpoint file
 amr.check_int       = 1000       # number of timesteps between checkpoints
@@ -48,12 +55,11 @@ amr.check_int       = 1000       # number of timesteps between checkpoints
 # PLOTFILES
 erf.plot_file_1         = plt        # root name of plotfile
 erf.plot_int_1          = 500         # number of timesteps between plotfiles
-erf.plot_vars_1         = density x_velocity y_velocity z_velocity pressure theta temp vorticity_x vorticity_y vorticity_z qv qc qrain rain_accum 
+erf.plot_vars_1         = density x_velocity y_velocity z_velocity pressure theta temp vorticity_x vorticity_y vorticity_z qv qc qrain qsat rain_accum 
 
 # SOLVER CHOICE
 erf.use_gravity = true
 erf.buoyancy_type = 1
-erf.use_coriolis = false
 
 #erf.les_type = "Smagorinsky"
 erf.Cs              = 0.25
@@ -68,16 +74,3 @@ erf.alpha_C           = 200.33
 erf.moisture_model = "Kessler"
 erf.use_moist_background = true
 
-# PROBLEM PARAMETERS (optional)
-prob.z_tr = 12000.0
-prob.height = 1200.0
-prob.theta_0 = 300.0
-prob.theta_tr = 343.0
-prob.T_tr = 213.0
-prob.x_c = 0.0
-prob.y_c = 0.0
-prob.z_c = 2000.0
-prob.x_r = 10000.0
-prob.y_r = 10000.0
-prob.z_r = 2000.0
-prob.theta_c = 0.0

Exec/DevTests/Hurricane/inputs_20210819_Henri → Exec/DevTests/Hurricane/inputs_20200826_Laura_NoAMR

@@ -1,5 +1,5 @@
 # ------------------  INPUTS TO MAIN PROGRAM  -------------------
-max_step = 51840
+max_step = 1000
 #max_step = 1
 stop_time = 1000000.0
 
@@ -8,11 +8,11 @@ amrex.fpe_trap_invalid = 1
 fabarray.mfiter_tile_size = 2048 1024 2048
 
 # PROBLEM SIZE & GEOMETRY
-geometry.prob_lo     = 1000001 2000001    0.
-geometry.prob_hi     = 4249999 4999999 25000.
-amr.n_cell           =  500    500    100    # dx=dy=dz=250 m
+geometry.prob_lo     = -1999999  1499999    0.
+geometry.prob_hi     =  3999999  5259999 25000.
+amr.n_cell           =  1000    500    100    # dx=dy=dz=250 m
 
-erf.IC_file = "ERF_IC_Henri_gdas1.fnl0p25.2021081906.f00.bin"
+erf.IC_file = "ERF_IC_Laura_LargeDomain.bin"
 
 # periodic in x to match WRF setup
 # - as an alternative, could use symmetry at x=0 and outflow at x=25600
@@ -24,6 +24,11 @@ yhi.type = "Open"
 zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
+erf.use_coriolis = true
+erf.coriolis_3d = true
+erf.latitude = 30.0
+erf.rotational_time_period = 86165.734375
+
 #erf.rayleigh_damp_T = true
 #erf.rayleigh_zdamp = 5000.0
 #erf.rayleigh_dampcoef = 0.005
@@ -47,13 +52,12 @@ amr.check_int       = 1000       # number of timesteps between checkpoints
 
 # PLOTFILES
 erf.plot_file_1         = plt        # root name of plotfile
-erf.plot_int_1          = 500         # number of timesteps between plotfiles
+erf.plot_int_1          = 1000         # number of timesteps between plotfiles
 erf.plot_vars_1         = density x_velocity y_velocity z_velocity pressure theta temp vorticity_x vorticity_y vorticity_z qv qc qrain rain_accum 
 
 # SOLVER CHOICE
 erf.use_gravity = true
 erf.buoyancy_type = 1
-erf.use_coriolis = false
 
 #erf.les_type = "Smagorinsky"
 erf.Cs              = 0.25
@@ -68,16 +72,3 @@ erf.alpha_C           = 200.33
 erf.moisture_model = "Kessler"
 erf.use_moist_background = true
 
-# PROBLEM PARAMETERS (optional)
-prob.z_tr = 12000.0
-prob.height = 1200.0
-prob.theta_0 = 300.0
-prob.theta_tr = 343.0
-prob.T_tr = 213.0
-prob.x_c = 0.0
-prob.y_c = 0.0
-prob.z_c = 2000.0
-prob.x_r = 10000.0
-prob.y_r = 10000.0
-prob.z_r = 2000.0
-prob.theta_c = 0.0

Exec/DevTests/Hurricane/inputs_20200826_Laura → Exec/DevTests/Hurricane/inputs_20210819_Henri_NoAMR

@@ -8,11 +8,11 @@ amrex.fpe_trap_invalid = 1
 fabarray.mfiter_tile_size = 2048 1024 2048
 
 # PROBLEM SIZE & GEOMETRY
-geometry.prob_lo     = -749999 1249999    0.
-geometry.prob_hi     = 2249999 4999999 25000.
-amr.n_cell           =  500    500    100    # dx=dy=dz=250 m
+geometry.prob_lo     = -1299999 2500001    0.
+geometry.prob_hi     =  5999999 5999999 25000.
+amr.n_cell           =  1000    500    100    # dx=dy=dz=250 m
 
-erf.IC_file = "ERF_IC_Laura_2020082600.bin"
+erf.IC_file = "ERF_IC_HenriERA5_20210819_VeryLarge.bin"
 
 # periodic in x to match WRF setup
 # - as an alternative, could use symmetry at x=0 and outflow at x=25600
@@ -24,6 +24,11 @@ yhi.type = "Open"
 zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
+erf.use_coriolis = true
+erf.coriolis_3d = true
+erf.latitude = 30.0
+erf.rotational_time_period = 86165.734375
+
 #erf.rayleigh_damp_T = true
 #erf.rayleigh_zdamp = 5000.0
 #erf.rayleigh_dampcoef = 0.005
@@ -43,7 +48,7 @@ amr.max_level       = 0       # maximum level number allowed
 # CHECKPOINT FILES
 amr.check_file      = chk        # root name of checkpoint file
 amr.check_int       = 1000       # number of timesteps between checkpoints
-#amr.restart         = chk15000
+#amr.restart         = chk47000
 
 # PLOTFILES
 erf.plot_file_1         = plt        # root name of plotfile
@@ -53,7 +58,6 @@ erf.plot_vars_1         = density x_velocity y_velocity z_velocity pressure thet
 # SOLVER CHOICE
 erf.use_gravity = true
 erf.buoyancy_type = 1
-erf.use_coriolis = false
 
 #erf.les_type = "Smagorinsky"
 erf.Cs              = 0.25
@@ -65,19 +69,5 @@ erf.dynamic_viscosity = 200.33 # [kg/(m-s)]
 erf.alpha_T           = 200.33 # [m^2/s]
 erf.alpha_C           = 200.33
 
-erf.moisture_model = "Kessler"
+erf.moisture_model = "None"
 erf.use_moist_background = true
-
-# PROBLEM PARAMETERS (optional)
-prob.z_tr = 12000.0
-prob.height = 1200.0
-prob.theta_0 = 300.0
-prob.theta_tr = 343.0
-prob.T_tr = 213.0
-prob.x_c = 0.0
-prob.y_c = 0.0
-prob.z_c = 2000.0
-prob.x_r = 10000.0
-prob.y_r = 10000.0
-prob.z_r = 2000.0
-prob.theta_c = 0.0