fix anelastic and add heffte as option (still a WIP) (#1838)

* fix anelastic and add heffte as option (still a WIP)

* fix for gpu

Exec/Make.ERF

@@ -181,10 +181,20 @@ VPATH_LOCATIONS += $(HEFFTE_HOME)/include
 INCLUDE_LOCATIONS += $(HEFFTE_HOME)/include
 LIBRARY_LOCATIONS += $(HEFFTE_HOME)/lib
 LIBRARIES += -lheffte
+ifeq ($(USE_CUDA),TRUE)
+ libraries += -lcufft
+else ifeq ($(USE_HIP),TRUE)
+ # Use rocFFT. ROC_PATH is defined in amrex
+ INCLUDE_LOCATIONS += $(ROC_PATH)/rocfft/include
+ LIBRARY_LOCATIONS += $(ROC_PATH)/rocfft/lib
+ LIBRARIES += -L$(ROC_PATH)/rocfft/lib -lrocfft
+else
+ libraries += -lfftw3_mpi -lfftw3f -lfftw3
+endif
 endif
 
 ifeq ($(USE_WW3_COUPLING), TRUE)
- DEFINES += -DERF_USE_WW3_COUPLING
+DEFINES += -DERF_USE_WW3_COUPLING
 endif
 
 ERF_LSM_DIR = $(ERF_SOURCE_DIR)/LandSurfaceModel

Exec/RegTests/DensityCurrent/inputs_anelastic

@@ -1,12 +1,8 @@
 # ------------------ INPUTS TO MAIN PROGRAM -------------------
-max_step = 999999
 stop_time = 900.0
 
 erf.anelastic = 1
-erf.no_substepping = 1
-erf.buoyancy_type = 2
 
-erf.use_terrain = true
 erf.use_terrain = false
 
 amrex.fpe_trap_invalid = 1
@@ -15,32 +11,27 @@ fabarray.mfiter_tile_size = 1024 1024 1024
 
 # PROBLEM SIZE & GEOMETRY
 
-#geometry.prob_lo = -12800. 0. 0.
-#geometry.prob_hi = 12800. 1600. 6400.
-#xlo.type = "Outflow"
-
 geometry.prob_lo = 0. 0. 0.
-geometry.prob_hi = 25600. 1600. 6400.
-xlo.type = "Symmetry"
-
-xhi.type = "HO_Outflow"
+geometry.prob_hi = 25600. 100. 6400.
 
 geometry.is_periodic = 0 1 0
 
-amr.n_cell = 256 16 64 # dx=dy=dz=100 m, Straka et al 1993 / Xue et al 2000
+amr.n_cell = 256 1 64 # dx=dy=dz=100 m, Straka et al 1993 / Xue et al 2000
+
+xlo.type = "Symmetry"
+xhi.type = "HO_Outflow"
 
 zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
 erf.fixed_dt = 1.0 # fixed time step [s] -- Straka et al 1993
-erf.fixed_fast_dt = 0.25 # fixed time step [s] -- Straka et al 1993
 
 # DIAGNOSTICS & VERBOSITY
 erf.sum_interval = 1 # timesteps between computing mass
-erf.v = 0 # verbosity in ERF.cpp
-erf.mg_v = 1 # verbosity in ERF.cpp
+erf.v = 1 # verbosity in ERF.cpp
 amr.v = 1 # verbosity in Amr.cpp
+erf.mg_v = 1 # verbosity in ERF.cpp
 
 # REFINEMENT / REGRIDDING
 amr.max_level = 0 # maximum level number allowed

Source/DataStructs/ERF_DataStruct.H

@@ -139,9 +139,6 @@ struct SolverChoice {
 
  // Which expression (1,2 or 3) to use for buoyancy
  pp.query("buoyancy_type", buoyancy_type);
- if (buoyancy_type != 1 && buoyancy_type != 2 && buoyancy_type != 3 && buoyancy_type != 4) {
- amrex::Abort("buoyancy_type must be 1, 2, 3 or 4");
- }
 
  // What type of land surface model to use
  static std::string lsm_model_string = "None";
@@ -168,17 +165,14 @@ struct SolverChoice {
  // Use lagged_delta_rt in the fast integrator?
  pp.query("use_lagged_delta_rt", use_lagged_delta_rt);
 
- if (!use_lagged_delta_rt && !(terrain_type == TerrainType::Moving)) {
- amrex::Error("Can't turn off lagged_delta_rt when terrain not moving");
- }
-
  // These default to true but are used for unit testing
  pp.query("use_gravity", use_gravity);
  gravity = use_gravity? CONST_GRAV: 0.0;
 
  pp.query("c_p", c_p);
  rdOcp = R_d / c_p;
 
+ read_int_string(max_level, "no_substepping", no_substepping, 0);
 
 #if defined(ERF_USE_POISSON_SOLVE)
  // Should we project the initial velocity field to make it divergence-free?
@@ -194,6 +188,7 @@ struct SolverChoice {
  // If anelastic at all, we do not advect rho -- it is always == rho0
  if (any_anelastic == 1) {
  constant_density = true;
+ buoyancy_type = 2; // uses Tprime
  } else {
  constant_density = false; // We default to false but allow the user to set it
  pp.query("constant_density", constant_density);
@@ -203,24 +198,20 @@ struct SolverChoice {
  pp.query("ncorr", ncorr);
  pp.query("poisson_abstol", poisson_abstol);
  pp.query("poisson_reltol", poisson_reltol);
-#else
- anelastic.resize(max_level+1);
- for (int i = 0; i <= max_level; ++i) anelastic[i] = 0;
-#endif
 
- read_int_string(max_level, "no_substepping", no_substepping, 0);
-
- pp.query("force_stage1_single_substep", force_stage1_single_substep);
-
-#if defined(ERF_USE_POISSON_SOLVE)
  for (int lev = 0; lev <= max_level; lev++) {
- if (anelastic[lev] != 0 && no_substepping[lev] == 0)
+ if (anelastic[lev] != 0)
  {
  no_substepping[lev] = 1;
  }
  }
+#else
+ anelastic.resize(max_level+1);
+ for (int i = 0; i <= max_level; ++i) anelastic[i] = 0;
 #endif
 
+ pp.query("force_stage1_single_substep", force_stage1_single_substep);
+
  // Include Coriolis forcing?
  pp.query("use_coriolis", use_coriolis);
 
@@ -344,12 +335,6 @@ struct SolverChoice {
  }
  }
 
- // Warn for PBL models and moisture - these may not yet be compatible
- for (int lev = 0; lev <= max_level; lev++) {
- if ((moisture_type != MoistureType::None) && (turbChoice[lev].pbl_type != PBLType::None)) {
- amrex::Warning("\n*** WARNING: Moisture may not yet be compatible with PBL models, \n proceed with caution ***");
- }
- }
 
  // Which type of refinement
  static std::string coupling_type_string = "TwoWay";
@@ -401,30 +386,56 @@ struct SolverChoice {
  // turbine. ie. the sampling distance will be this number multiplied
  // by the diameter of the turbine
  pp.query("sampling_distance_by_D", sampling_distance_by_D);
- if(windfarm_type==WindFarmType::SimpleAD and sampling_distance_by_D < 0.0) {
+ pp.query("turb_disk_angle_from_x", turb_disk_angle);
+
+ pp.query("windfarm_x_shift",windfarm_x_shift);
+ pp.query("windfarm_y_shift",windfarm_y_shift);
+ // Test if time averaged data is to be output
+ pp.query("time_avg_vel",time_avg_vel);
+
+ check_params(max_level);
+ }
+
+ void check_params(int max_level)
+ {
+ // Warn for PBL models and moisture - these may not yet be compatible
+ for (int lev = 0; lev <= max_level; lev++) {
+ if ((moisture_type != MoistureType::None) && (turbChoice[lev].pbl_type != PBLType::None)) {
+ amrex::Warning("\n*** WARNING: Moisture may not yet be compatible with PBL models, \n proceed with caution ***");
+ }
+ }
+ //
+ // Buoyancy type check
+ //
+ if (buoyancy_type != 1 && buoyancy_type != 2 && buoyancy_type != 3 && buoyancy_type != 4) {
+ amrex::Abort("buoyancy_type must be 1, 2, 3 or 4");
+ }
+
+ if (!use_lagged_delta_rt && !(terrain_type == TerrainType::Moving)) {
+ amrex::Error("Can't turn off lagged_delta_rt when terrain not moving");
+ }
+
+ //
+ // Wind farm checks
+ //
+ if (windfarm_type==WindFarmType::SimpleAD and sampling_distance_by_D < 0.0) {
  amrex::Abort("To use simplified actuator disks, you need to provide a variable"
  " erf.sampling_distance_by_D in the inputs which specifies the upstream"
  " distance as a factor of the turbine diameter at which the incoming free stream"
  " velocity will be computed at.");
  }
- pp.query("turb_disk_angle_from_x", turb_disk_angle);
- if(windfarm_type==WindFarmType::SimpleAD and turb_disk_angle < 0.0) {
+ if (windfarm_type==WindFarmType::SimpleAD and turb_disk_angle < 0.0) {
  amrex::Abort("To use simplified actuator disks, you need to provide a variable"
  " erf.turb_disk_angle_from_x in the inputs which is the angle of the face of the"
  " turbine disk from the x-axis. A turbine facing an oncoming flow in the x-direction"
  " will have turb_disk_angle value of 90 deg.");
  }
-
- pp.query("windfarm_x_shift",windfarm_x_shift);
- pp.query("windfarm_y_shift",windfarm_y_shift);
- if(windfarm_loc_type == WindFarmLocType::lat_lon and (windfarm_x_shift < 0.0 or windfarm_y_shift < 0.0)) {
+ if (windfarm_loc_type == WindFarmLocType::lat_lon and (windfarm_x_shift < 0.0 or windfarm_y_shift < 0.0)) {
  amrex::Abort("You are using windfarms with latitude-logitude option to position the turbines."
  " For this you should provide the inputs erf.windfarm_x_shift and"
  " erf.windfarm_y_shift which are the values by which the bounding box of the"
  " windfarm is shifted from the x and the y axes.");
  }
- // Test if time averaged data is to be output
- pp.query("time_avg_vel",time_avg_vel);
  }
 
  void display(int max_level)

Source/ERF.H

@@ -137,6 +137,8 @@ public:
 #ifdef ERF_USE_POISSON_SOLVE
  // Project the velocities to be divergence-free
  void project_velocities (int lev, amrex::Real dt, amrex::Vector<amrex::MultiFab >& vars, amrex::MultiFab& p);
+ void solve_with_heffte (int lev, amrex::MultiFab& rhs, amrex::MultiFab& soln,
+ amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>& fluxes);
 
  // Project the velocities to be divergence-free with a thin body
  void project_velocities_tb (int lev, amrex::Real dt, amrex::Vector<amrex::MultiFab >& vars, amrex::MultiFab& p);
@@ -940,6 +942,7 @@ private:
  static int verbose;
 #ifdef ERF_USE_POISSON_SOLVE
  static int mg_verbose;
+ static bool use_heffte;
 #endif
 
  // Diagnostic output interval

Source/ERF.cpp

@@ -41,7 +41,8 @@ int ERF::fixed_mri_dt_ratio = 0;
 // Dictate verbosity in screen output
 int ERF::verbose = 0;
 #ifdef ERF_USE_POISSON_SOLVE
-int ERF::mg_verbose = 0;
+int ERF::mg_verbose = 0;
+bool ERF::use_heffte = false;
 #endif
 
 // Frequency of diagnostic output
@@ -1380,6 +1381,7 @@ ERF::ReadParameters ()
  pp.query("v", verbose);
 #ifdef ERF_USE_POISSON_SOLVE
  pp.query("mg_v", mg_verbose);
+ pp.query("use_heffte", use_heffte);
 #endif
 
  // Frequency of diagnostic output

Source/SourceTerms/ERF_Src_headers.H

@@ -21,7 +21,9 @@ void make_buoyancy (amrex::Vector< amrex::MultiFab>& S_data,
  const amrex::Geometry geom,
  const SolverChoice& solverChoice,
  const amrex::MultiFab* r0,
- const int n_qstate);
+ const amrex::MultiFab* p0,
+ const int n_qstate,
+ const int anelastic);
 
 void make_sources (int level, int nrk,
  amrex::Real dt, amrex::Real time,
@@ -55,6 +57,7 @@ void make_mom_sources (int level, int nrk,
  amrex::MultiFab& ymom_source,
  amrex::MultiFab& zmom_source,
  amrex::MultiFab* r0,
+ amrex::MultiFab* p0,
  const amrex::Geometry geom,
  const SolverChoice& solverChoice,
  std::unique_ptr<amrex::MultiFab>& mapfac_m,