diff --git a/Exec/DevTests/Hurricane/ERF_Prob.cpp b/Exec/DevTests/Hurricane/ERF_Prob.cpp index 361346983..0e0f09cdc 100644 --- a/Exec/DevTests/Hurricane/ERF_Prob.cpp +++ b/Exec/DevTests/Hurricane/ERF_Prob.cpp @@ -1,6 +1,8 @@ #include "ERF_Prob.H" #include +#include "ERF_Interpolation_Bilinear.H" + using namespace amrex; std::unique_ptr @@ -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 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, diff --git a/Exec/DevTests/Hurricane/README.md b/Exec/DevTests/Hurricane/README.md index 023975e59..51aecd6ba 100644 --- a/Exec/DevTests/Hurricane/README.md +++ b/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`. diff --git a/Exec/DevTests/Hurricane/inputs_20210821_Henri b/Exec/DevTests/Hurricane/inputs_20200826_Laura_2LevelAMR similarity index 72% rename from Exec/DevTests/Hurricane/inputs_20210821_Henri rename to Exec/DevTests/Hurricane/inputs_20200826_Laura_2LevelAMR index bcc3c12e5..70b0cee41 100644 --- a/Exec/DevTests/Hurricane/inputs_20210821_Henri +++ b/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 diff --git a/Exec/DevTests/Hurricane/inputs_20210819_Henri b/Exec/DevTests/Hurricane/inputs_20200826_Laura_NoAMR similarity index 74% rename from Exec/DevTests/Hurricane/inputs_20210819_Henri rename to Exec/DevTests/Hurricane/inputs_20200826_Laura_NoAMR index a993ed080..c96730547 100644 --- a/Exec/DevTests/Hurricane/inputs_20210819_Henri +++ b/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 diff --git a/Exec/DevTests/Hurricane/inputs_20200826_Laura b/Exec/DevTests/Hurricane/inputs_20210819_Henri_NoAMR similarity index 76% rename from Exec/DevTests/Hurricane/inputs_20200826_Laura rename to Exec/DevTests/Hurricane/inputs_20210819_Henri_NoAMR index 859268fe4..e4c33e550 100644 --- a/Exec/DevTests/Hurricane/inputs_20200826_Laura +++ b/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 diff --git a/Source/Utils/ERF_Interpolation_Bilinear.H b/Source/Utils/ERF_Interpolation_Bilinear.H new file mode 100644 index 000000000..2b3d77f69 --- /dev/null +++ b/Source/Utils/ERF_Interpolation_Bilinear.H @@ -0,0 +1,74 @@ +#ifndef ERF_INTERPOLATE_BILINEAR_H_ +#define ERF_INTERPOLATE_BILINEAR_H_ + +#include "ERF_DataStruct.H" + +/** + * Bilinear interpolation routines. There are different versions for + * different scenarios that take in a different argument list + */ + +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 amrex::Real* xvec, const amrex::Real* yvec, const amrex::Real* zvec, + const amrex::Real dxvec, const amrex::Real dyvec, + const int nx, const int ny, const int nz, + const amrex::Real x, const amrex::Real y, const amrex::Real z, + const amrex::Real* varvec, + amrex::Real& tmp_var) +{ + int iloc=-1, jloc=-1, kloc=-1; + for(int k=0;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){ + } + + amrex::Real xlo = xvec[0] + iloc*dxvec; + amrex::Real ylo = yvec[0] + jloc*dyvec; + amrex::Real zlo = zvec[kloc]; + + amrex::Real w_x = (x - xlo)/dxvec; + amrex::Real w_y = (y - ylo)/dyvec; + amrex::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"; +} + +#endif diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index ab50e180a..c3c07dc99 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -5,6 +5,7 @@ CEXE_headers += ERF_Interpolation_WENO.H CEXE_headers += ERF_Interpolation_WENO_Z.H CEXE_headers += ERF_Interpolation.H CEXE_headers += ERF_Interpolation_1D.H +CEXE_headers += ERF_Interpolation_Bilinear.H CEXE_headers += ERF_Microphysics_Utils.H CEXE_headers += ERF_TerrainMetrics.H CEXE_headers += ERF_TileNoZ.H