diff --git a/Docs/sphinx_doc/MOST.rst b/Docs/sphinx_doc/MOST.rst index e644529a1..cbd33578d 100644 --- a/Docs/sphinx_doc/MOST.rst +++ b/Docs/sphinx_doc/MOST.rst @@ -149,13 +149,13 @@ In ERF, when the MOST boundary condition is applied, velocity and temperature in MOST Inputs ~~~~~~~~~~~~~~~~~~~ -To evaluate the fluxes with MOST, the surface rougness parameter :math:`z_{0}` must be specified. This quantity may be considered a constant or may be parameterized through the friction velocity :math:`u_{\star}`. ERF supports three methods for parameterizing the surface roughness: ``constant``, ``charnock``, and ``modified_charnock``. The latter two methods parameterize :math:`z_{0} = f(u_{\star})` and are described in `Jimenez & Dudhia, American Meteorological Society, 2018 `_. The rougness calculation method may be specified with +To evaluate the fluxes with MOST, the surface rougness parameter :math:`z_{0}` must be specified. This quantity may be considered a constant or may be parameterized through the friction velocity :math:`u_{\star}`. ERF supports four methods for parameterizing the surface roughness: ``constant``, ``charnock``, ``modified_charnock``, and ``wave_coupled``. The latter three methods parameterize :math:`z_{0} = f(u_{\star})` and are described in `Jimenez & Dudhia, American Meteorological Society, 2018 `_ and `Warner et. al, Ocean Modelling, 2010 `_. The rougness calculation method may be specified with :: - erf.most.roughness_type = STRING #Z_0 type (constant, charnock, modified_charnock) + erf.most.roughness_type = STRING #Z_0 type (constant, charnock, modified_charnock, wave_couples) -If the ``charnock`` method is employed, the :math:`a` constant may be specified with ``erf.most.charnock_constant`` (defaults to 0.0185). If the ``modified_charnock`` method is employed, the depth :math:`d` may be specified with ``erf.most.modified_charnock_depth`` (defaults to 30 m). +If the ``charnock`` method is employed, the :math:`a` constant may be specified with ``erf.most.charnock_constant`` (defaults to 0.0185). If the ``modified_charnock`` method is employed, the depth :math:`d` may be specified with ``erf.most.modified_charnock_depth`` (defaults to 30 m). If the ``wave_coupled`` method is employed, the user must provide wave height and mean wavelength data. When computing an average :math:`\overline{\phi}` for the MOST boundary, where :math:`\phi` denotes a generic variable, ERF supports a variety of approaches. Specifically, ``planar averages`` and ``local region averages`` may be computed with or without ``time averaging``. With each averaging methodology, the query point :math:`z` may be determined from the following procedures: specified vertical distance :math:`z_{ref}` from the bottom surface, specified :math:`k_{index}`, or (when employing terrain-fit coordinates) specified normal vector length :math:`z_{ref}`. The available inputs to the MOST boundary and their associated data types are diff --git a/Docs/sphinx_doc/building.rst b/Docs/sphinx_doc/building.rst index 481d8d4cd..1c709cb21 100644 --- a/Docs/sphinx_doc/building.rst +++ b/Docs/sphinx_doc/building.rst @@ -263,7 +263,7 @@ For Perlmutter at NERSC, look at the general instructions for building ERF using module load PrgEnv-gnu module load cudatoolkit -Then build ERF as, for example (specify your own path to the AMReX submodule in `ERF/Submodules/AMReX`): +Then build ERF as, for example (specify your own path to the AMReX submodule in ``ERF/Submodules/AMReX``): :: @@ -308,5 +308,47 @@ Finally, you can prepare your SLURM job script, using the following as a guide: ./ERF3d.gnu.MPI.CUDA.ex inputs_wrf_baseline max_step=100 ${GPU_AWARE_MPI}" \ > test.out -To submit your job script, do `sbatch [your job script]` and you can check its status by doing `squeue -u [your username]`. +To submit your job script, do ``sbatch [your job script]`` and you can check its status by doing ``squeue -u [your username]``. + +Kestrel (NREL) +~~~~~~~~~~~~~~ + +The `Kestrel `_ cluster is an HPE Cray machine +composed primarily of CPU compute nodes with 104 core +Intel Xeon Sapphire Rapids nodes. It also contains a GPU partition with 4 Nvidia H100 GPUs per node. + +As with Perlmutter, the GNU Make build system is preferred. To compile and run on CPUs, the default modules +loaded when logging into Kestrel can be used. If you are unsure about your environment, you can reset to +the default modules: :: + + module restore + +Then, build ERF using the cray compilers (if wishing to use other compilers, you can swap the ``PrgEnv-cray`` module +for another module as appropriate, see Kestrel user documentation for more details): :: + + make realclean; make -j COMP=cray + +For compiling and running on GPUs, the following commands can be used to set up your environment: :: + + module restore; module load PrgEnv-gnu/8.5.0; module load cray-libsci/23.05.1.4; module load cmake; module load cuda/12.3; module load cray-mpich/8.1.28; module load craype/2.7.30; + +And then compile: :: + + make realclean; make -j COMP=gnu USE_CUDA=TRUE + +When running on Kestrel, GPU node hours are charged allocation units (AUs) at 10 times the rate of CPU node hours. +For ERF, the performance running on a Kestrel GPU node with 4 GPUs is typically 10-20x running on a CPU node +with 96-104 MPI ranks per node, so the performance gain from on on GPUs is likely worth the higher charge +rate for node hours, in addition to providing faster time to solution. However, for smaller problem sizes, +or problems distributed across too many nodes (resulting in fewer than around 1 million cells/GPU), +the compute capability of the GPUs may be unsaturated and the performance gain from running on GPUs +may not justify the higher AU charge. The trade-off is problem dependent, so users may wish to assess +performance for their particular case and objectives in terms of wall time, AUs used, etc to determine the +optimal strategy if running large jobs. + +Another note about using Kestrel is that partial node allocations are possible, which means the full memory +available on each node may not be assigned by default. In general, using the ``--exclusive`` flag when +requesting nodes through the slurm scheduler, which will allocate entire nodes exlcusively for your request, +is recommended. Otherwise, memory intensive operations such as CUDA compilation may fail. You can alternatively +request a particular amount of memory with the ``--mem=XXX`` or ``--mem-per-cpu=XXX`` slurm inputs. diff --git a/Docs/sphinx_doc/figures/WindTurbine_EWP.png b/Docs/sphinx_doc/figures/WindTurbine_EWP.png new file mode 100644 index 000000000..15c006b4a Binary files /dev/null and b/Docs/sphinx_doc/figures/WindTurbine_EWP.png differ diff --git a/Docs/sphinx_doc/figures/WindTurbine_Fitch.png b/Docs/sphinx_doc/figures/WindTurbine_Fitch.png index 226ceb1cb..f66d661dd 100644 Binary files a/Docs/sphinx_doc/figures/WindTurbine_Fitch.png and b/Docs/sphinx_doc/figures/WindTurbine_Fitch.png differ diff --git a/Docs/sphinx_doc/theory/WindFarmModels.rst b/Docs/sphinx_doc/theory/WindFarmModels.rst index d9dc39182..ed709e697 100644 --- a/Docs/sphinx_doc/theory/WindFarmModels.rst +++ b/Docs/sphinx_doc/theory/WindFarmModels.rst @@ -4,7 +4,9 @@ Wind farm models Introduction ------------- -ERF supports models for wind farm parametrization in which the effects of wind turbines are represented by imposing a momentum sink on the mean flow and/or turbulent kinetic energy (TKE). Currently only the Fitch model (`Fitch et al. 2012`_) is supported. +ERF supports models for wind farm parametrization in which the effects of wind turbines are represented by imposing a momentum sink on the mean flow and/or turbulent kinetic energy (TKE). Currently only the Fitch model (`Fitch et al. 2012`_) and Explicit Wake Parametrization (EWP) model (`Volker et al. 2015`_) are supported. + +.. _Fitch model: Fitch model ------------ @@ -17,26 +19,27 @@ The Fitch model for wind farms introduced in `Fitch et al. 2012`_ models the ef \frac{\partial u_{ijk}}{\partial t} &= \frac{u_{ijk}}{|V|_{ijk}}\frac{\partial |V|_{ijk}}{\partial t} \\ \frac{\partial v_{ijk}}{\partial t} &= \frac{v_{ijk}}{|V|_{ijk}}\frac{\partial |V|_{ijk}}{\partial t} \\ - \frac{\partial \text{TKE}_{ijk}}{\partial t} &= \frac{0.5N^{ij}C_{TKE}(|V|_{ijk})|V|_{ijk}^3A_{ijk}}{z_{k+1}-z_k} + \frac{\partial \text{TKE}_{ijk}}{\partial t} &= \frac{0.5N_{ij}C_{TKE}(|V|_{ijk})|V|_{ijk}^3A_{ijk}}{\Delta x \Delta y (z_{k+1}-z_k)} where .. math:: - \frac{\partial |V|_{ijk}}{\partial t} = \frac{0.5N^{ij}C_T(|V|_{ijk})|V|_{ijk}^2A_{ijk}}{z_{k+1}-z_k} + \frac{\partial |V|_{ijk}}{\partial t} = \frac{0.5N_{ij}C_T(|V|_{ijk})|V|_{ijk}^2A_{ijk}}{\Delta x \Delta y (z_{k+1}-z_k)} + +where `u` and `v` are horizontal components of velocity, `|V|` is the velocity magnitude, :math:`N_{ij}` is the number of turbines in cell :math:`(i,j)` (see Fig. :numref:`fig:WindFarm`), :math:`C_T` is the coefficient of thrust of the turbines, :math:`C_{TKE}` is the fraction of energy converted to turbulent kinetic energy -- both of which are functions of the velocity magnitude, and :math:`A_{ijk}` is the area intersected by the swept area of the turbine between levels :math:`z=z_k` and :math:`z= z_{k+1}`, and is explained in the next section. -where `u` and `v` are horizontal components of velocity, `|V|` is the velocity magnitude, :math:`N^{ij}` is the number of turbines in cell :math:`(i,j)`, :math:`C_T` is the coefficient of thrust of the turbines, :math:`C_{TKE}` is the fraction of energy converted to turbulent kinetic energy -- both of which are functions of the velocity magnitude, and :math:`A_{ijk}` is the area intersected by the swept area of the turbine between levels :math:`z=z_k` and :math:`z= z_{k+1}`, and is explained in the next section. Intersected area :math:`A_{ijk}` -_________________________________ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Consider :math:`A_k^{k+1}` -- the area intersected by the swept area of the wind turbine between :math:`z=z_k` and :math:`z = z_{k+1}`. We have (see Figs. `2` and `3` below) +Consider :math:`A_k^{k+1}` -- the area intersected by the swept area of the wind turbine between :math:`z=z_k` and :math:`z = z_{k+1}`. We have (see Figs. :numref:`fig:WindTurbine_Fitch` and :numref:`fig:Fitch_Aijk` below) .. math:: A_k = \frac{\pi R^2}{2} - A_{ks} -where :math:`A_{ks}` is the area of the segment of the circle as shown in Fig. `3`. We have from geometry, :math:`d_k = \min(|z_k - z_c|,R)` is the perpendicular distance of the center of the turbine to :math:`z = z_k`, where :math:`z_c` is the height of the center of the turbine from the ground. The area of the segment is +where :math:`A_{ks}` is the area of the segment of the circle as shown in Fig. :numref:`fig:Fitch_Aijk`. We have from geometry, :math:`d_k = \min(|z_k - h|,R)` is the perpendicular distance of the center of the turbine to :math:`z = z_k`, where :math:`h` is the height of the center of the turbine from the ground. The area of the segment is .. math:: @@ -50,13 +53,13 @@ Hence, we have the intersected area :math:`A_{ijk}\equiv A_k^{k+1}` as A_k^{k+1} = \begin{cases} - |A_k - A_{k+1}| & \text{if } (z_k - z_c)(z_{k+1}-z_c) > 0 \\ - |A_k + A_{k+1}| & \text{if } (z_k - z_c)(z_{k+1}-z_c) \le 0 \\ + |A_k - A_{k+1}| & \text{if } (z_k - h)(z_{k+1}-h) > 0 \\ + |A_k + A_{k+1}| & \text{if } (z_k - h)(z_{k+1}-h) \le 0 \\ \end{cases} An example of the Fitch model is in ``Exec/Fitch`` -.. 1: +.. _fig:WindFarm: .. figure:: ../figures/WindFarm_Fitch.png :width: 300 @@ -64,18 +67,150 @@ An example of the Fitch model is in ``Exec/Fitch`` Horizontal view of the wind farm in the Fitch model -- shows a wind farm in cell `(i,j)` with 5 wind turbines. The turbines are discretized only in the vertical direction. -.. 2: +.. _fig:WindTurbine_Fitch: .. figure:: ../figures/WindTurbine_Fitch.png - :width: 300 + :width: 400 :align: center The vertical discretization of the wind turbine in the Fitch model. -.. 3: +.. _fig:Fitch_Aijk: .. figure:: ../figures/FitchModel_A_ijk.png :width: 400 :align: center The area terminology in the Fitch model. The circle represents the area swept by the turbine blades. + + +.. _explicit-wake-parametrization-ewp-model: + +Explicit Wake Parametrization (EWP) model +----------------------------------------- + +The Explicit Wake Parametrization (EWP) model [`Volker et al. 2015`_] is very similar to the Fitch model, and models the effect of wind farms as source terms in the governing equations for the horizontal components of momentum (i.e., :math:`x` and :math:`y` momentum) and the turbulent kinetic energy (TKE). At a given cell :math:`(i,j,k)`, the source terms in the governing equations are: + +.. math:: + \frac{\partial u_{ijk}}{\partial t} = -\sqrt{\frac{\pi}{8}}\frac{N_{ij}c_tr_0^2\overline{u}_0^2}{\Delta x \Delta y \sigma_e} + \exp\left\{-\frac{1}{2}\left(\frac{z-h}{\sigma_e}\right)^2\right\}\cos(\phi(k)) + +.. math:: + \frac{\partial v_{ijk}}{\partial t} = -\sqrt{\frac{\pi}{8}}\frac{N_{ij}c_tr_0^2\overline{u}_0^2}{\Delta x \Delta y \sigma_e} + \exp\left\{-\frac{1}{2}\left(\frac{z-h}{\sigma_e}\right)^2\right\}\sin(\phi(k)) + +.. math:: + \frac{\partial \text{TKE}_{ijk}}{\partial t} = -N_{ij}\rho A_rc_t\langle \overline{u}_{i,h}\overline{u'^2}_{i,h}\rangle + +with + +.. math:: + \sigma_e = \frac{\overline{u}_0}{3KL}\left[\left(\frac{2KL}{\overline{u}_0} + \sigma_0^2\right)^{\frac{3}{2}} - \sigma_0^3\right] + +where :math:`N_{ij}` is the number of turbines in cell :math:`(i,j)`, :math:`c_t` is the thrust coefficient, :math:`r_0` is the rotor radius, :math:`\overline{u}_0` is the mean advection velocity at hub height, :math:`h` is the hub height, :math:`\sigma_0 \approx 1.7 r_0` [`Volker et al. 2017`_] is a length scale that accounts for near-wake expansion, :math:`L` is the downstream distance that the wake travels within the cell approximated as a fraction of the cell size, :math:`K` is the turbulence eddy diffusivity (:math:`m^2/s`), :math:`\Delta x` and :math:`\Delta y` are the mesh sizes in the horizontal directions, and :math:`\phi(k)` is the wind direction with respect to the x-axis. :math:`\overline{u}_{i,h}` and :math:`\overline{u'}_{i,h}` are the mean and the fluctuating values of the velocity components (subscript :math:`i` is the direction index) at the hub height :math:`h`, :math:`A_r = \pi r^2` is the swept area of the rotor and :math:`\rho` is the density of air. + +The EWP model does not have a concept of intersected area by the turbine rotor like the Fitch model (see :ref:`Fitch model`). The exponential factor in the source terms for the velocities models the effect of the rotor for the momentum sinks (see Fig. :numref:`fig:WindTurbine_EWP`), and the turbulent kinetic energy source term only depends on the density, hub-height mean velocities and fluctuations, and the total swept area of the rotor :math:`A_r`. + +.. _fig:WindTurbine_EWP: + +.. figure:: ../figures/WindTurbine_EWP.png + :width: 400 + :align: center + + In the EWP model, the exponential factor in the source terms for the velocities models the effect of the rotor for the momentum sinks unlike the Fitch model which computes the + intersected area (see Fig. :numref:`fig:WindTurbine_Fitch`). + +.. _`Volker et al. 2015`: https://gmd.copernicus.org/articles/8/3715/2015/ +.. _`Volker et al. 2017`: https://iopscience.iop.org/article/10.1088/1748-9326/aa5d86 + + +.. _Inputs: + +Inputs for wind farm parametrization models +------------------------------------------------------------ + +Fitch, EWP +~~~~~~~~~~~ + +The following are the inputs required for simulations with Fitch, EWP models. + +.. code-block:: cpp + + // The parametrization model - Fitch, EWP + erf.windfarm_type = "Fitch" + + // How are the turbine locations specified? - using latitude-longitude + // format or x-y format? lat_lon or x_y + erf.windfarm_loc_type = "lat_lon" + + // If using lat_lon, then the latitude and longitude of + // the lower bottom corner of the domain has to be specified + // The following means 35 deg N, 100 deg W (note the negative sign) + erf.latitude_lo = 35.0 + erf.longitude_lo = -100.0 + + // Table having the wind turbine locations + erf.windfarm_loc_table = "windturbines_1WT.txt" + + // Table having the specifications of the wind turbines. All turbines are assumed to + // have the same specifications + erf.windfarm_spec_table = "wind-turbine_1WT.tbl" + +1. ``erf.windfarm_type`` has to be one of the supported models - ``Fitch``, ``EWP``. +2. ``erf.windfarm_loc_type`` is a variable to specify how the wind turbine locations in the wind farm is specified. If using the latitude and longitude of the turbine location, this has to be ``lat_lon`` or if using x and y co-ordinates to specify the turbine locations, this input is ``x_y``. + + - If using ``lat_lon`` format, ``erf.latitude_lo`` and ``erf.longitude_lo`` specifies the latitude and longitude of the lower bottom corner of the domain box. ie. if the domain box is specified as + + .. code-block:: cpp + + geometry.prob_lo = -25000. 0. -10000 + geometry.prob_hi = 25000. 10000. 10000. + + then ``(erf.latitude_lo, erf.longitude_lo)`` corresponds to ``(-25000, 0, -10000)``. + + - If using ``x_y`` format, there is no need to specify the ``erf.latitude_lo`` and ``erf.longitude_lo``. + +3. The ``erf.windfarm_loc_table`` specifies the locations of the wind turbines in the wind farm. + + - For the latitude-longitude format, an example is as below. Each line specifies the latitude and longitude of the wind turbine location. The third entry simply has to be always 1 (WRF requires the third entry to be always 1, so maintaining same format here). The first entry means that the turbine is located at ``35.7828 deg N, 99.0168 deg W`` (note the negative sign in the entry corresponding to West). + + .. code-block:: cpp + + 35.7828828829 -99.0168 1 + 35.8219219219 -99.0168 1 + 35.860960961 -99.0168 1 + 35.9 -99.0168 1 + 35.939039039 -99.0168 1 + 35.9780780781 -99.0168 1 + 36.0171171171 -99.0168 1 + 35.7828828829 -98.9705333333 1 + + - For the x-y format, an example is as below. Each line specifies the x and y co-ordinates of the wind turbine location in metres + + .. code-block:: cpp + + 89264.99080053 91233.3333309002 + 89259.1966417755 95566.6666710007 + 89254.1277665419 99900.0000000001 + 89249.7842982733 104233.333329 + 89246.1663427532 108566.6666691 + 89243.2739881117 112899.9999981 + 93458.6633652711 86900.0000019001 + 93450.4377452458 91233.3333309002 + 93442.9032518779 95566.6666710007 + +4. The ``erf.windfarm_spec_table`` gives the specifications of the wind turbines in the wind farm. All wind turbines are assumed to have the same specifications. Here is a sample specifications table. + +.. code-block:: cpp + + 4 + 119.0 178.0 0.130 2.0 + 9 0.805 50.0 + 10 0.805 50.0 + 11 0.805 50.0 + 12 0.805 50.0 + +The first line is the number of pairs of entries for the power curve and thrust coefficient (there are 4 entries in this table which are in the last four lines of the table). +The second line gives the height in meters of the turbine hub, the diameter in +meters of the rotor, the standing thrust coefficient, and the nominal power of the turbine in MW. +The remaining lines (four in this case) contain the three values of: wind speed (m/s), thrust coefficient, and power production in kW. diff --git a/Exec/DevTests/Bomex/input_Kessler b/Exec/DevTests/Bomex/input_Kessler index de31b3a20..422d9caf6 100644 --- a/Exec/DevTests/Bomex/input_Kessler +++ b/Exec/DevTests/Bomex/input_Kessler @@ -1,12 +1,12 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -stop_time = 21600 +stop_time = 43200 amrex.fpe_trap_invalid = 1 fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY -geometry.prob_extent = 6400 6400 4000 +geometry.prob_extent = 3200 3200 4000 amr.n_cell = 32 32 100 geometry.is_periodic = 1 1 0 @@ -17,20 +17,22 @@ erf.most.flux_type = "custom" erf.most.ustar = 0.28 # ustar erf.most.tstar = 8.0e-3 # theta flux erf.most.qstar = 5.2e-5 # qv flux -erf.most.z0 = 0.1 +erf.most.z0 = 0.1 erf.most.zref = 20.0 -zhi.type = "HO_Outflow" - +# NOTE: This should have a qv grad too (use hoextrapcc?!) +zhi.type = "SlipWall" +#zhi.theta_grad = 0.00365 + # TIME STEP CONTROL -erf.fixed_dt = 2.0 # fixed time step depending on grid resolution +erf.fixed_dt = 0.5 # fixed time step depending on grid resolution # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass erf.v = 1 # verbosity in ERF.cpp amr.v = 1 # verbosity in Amr.cpp erf.data_log = "surf" "mean" "flux" "subgrid" -erf.profile_int = 30 +erf.profile_int = 120 # REFINEMENT / REGRIDDING amr.max_level = 0 # maximum level number allowed @@ -41,16 +43,16 @@ erf.check_int = 600 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 30 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhotheta x_velocity y_velocity z_velocity pressure temp theta qt qp qv qc qi qrain +erf.plot_int_1 = 120 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta x_velocity y_velocity z_velocity pressure temp theta qt qp qv qc qi # SOLVER CHOICE erf.alpha_T = 0.0 -erf.alpha_C = 1.0 +erf.alpha_C = 0.0 erf.use_gravity = true -erf.dycore_horiz_adv_type = Upwind_5th -erf.dycore_vert_adv_type = Upwind_5th +erf.dycore_horiz_adv_type = Upwind_3rd +erf.dycore_vert_adv_type = Upwind_3rd erf.dryscal_horiz_adv_type = WENOZ5 erf.dryscal_vert_adv_type = WENOZ5 erf.moistscal_horiz_adv_type = WENOZ5 @@ -59,17 +61,25 @@ erf.moistscal_vert_adv_type = WENOZ5 erf.moisture_model = "Kessler" erf.molec_diff_type = "None" + erf.les_type = "Smagorinsky" -erf.Cs = 0.15 -#erf.Sc_t = 0.3333333 +erf.Cs = 0.17 + +#erf.les_type = "Deardorff" +#erf.Ck = 0.1 +#erf.sigma_k = 1.0 +#erf.Ce = 0.1 + +erf.Pr_t = 0.333333 +erf.Sc_t = 0.333333 erf.init_type = "input_sounding" erf.init_sounding_ideal = true -erf.add_custom_rhotheta_forcing = true -erf.add_custom_moisture_forcing = true -erf.add_custom_geostrophic_profile = true -erf.add_custom_w_subsidence = true +erf.add_custom_rhotheta_forcing = true +erf.add_custom_moisture_forcing = true +erf.add_custom_geostrophic_profile = true +erf.add_custom_w_subsidence = true erf.custom_forcing_uses_primitive_vars = true # Higher values of perturbations lead to instability diff --git a/Exec/DevTests/Bomex/input_SAM b/Exec/DevTests/Bomex/input_SAM index a108e3de5..5e70cf9fb 100644 --- a/Exec/DevTests/Bomex/input_SAM +++ b/Exec/DevTests/Bomex/input_SAM @@ -1,13 +1,13 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -stop_time = 21600 +stop_time = 43200 amrex.fpe_trap_invalid = 1 fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY -geometry.prob_extent = 6400 6400 4000 -amr.n_cell = 64 64 100 +geometry.prob_extent = 3200 3200 4000 +amr.n_cell = 32 32 100 geometry.is_periodic = 1 1 0 @@ -17,20 +17,22 @@ erf.most.flux_type = "custom" erf.most.ustar = 0.28 # ustar erf.most.tstar = 8.0e-3 # theta flux erf.most.qstar = 5.2e-5 # qv flux -erf.most.z0 = 0.1 +erf.most.z0 = 0.1 erf.most.zref = 20.0 +# NOTE: This should have a qv grad too (use hoextrapcc?!) zhi.type = "SlipWall" - +#zhi.theta_grad = 0.00365 + # TIME STEP CONTROL -erf.fixed_dt = 1.0 # fixed time step depending on grid resolution +erf.fixed_dt = 0.5 # fixed time step depending on grid resolution # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass erf.v = 1 # verbosity in ERF.cpp amr.v = 1 # verbosity in Amr.cpp erf.data_log = "surf" "mean" "flux" "subgrid" -erf.profile_int = 60 +erf.profile_int = 120 # REFINEMENT / REGRIDDING amr.max_level = 0 # maximum level number allowed @@ -41,16 +43,16 @@ erf.check_int = 600 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 60 # number of timesteps between plotfiles +erf.plot_int_1 = 120 # number of timesteps between plotfiles erf.plot_vars_1 = density rhotheta x_velocity y_velocity z_velocity pressure temp theta qt qp qv qc qi # SOLVER CHOICE erf.alpha_T = 0.0 -erf.alpha_C = 1.0 +erf.alpha_C = 0.0 erf.use_gravity = true -erf.dycore_horiz_adv_type = Upwind_5th -erf.dycore_vert_adv_type = Upwind_5th +erf.dycore_horiz_adv_type = Upwind_3rd +erf.dycore_vert_adv_type = Upwind_3rd erf.dryscal_horiz_adv_type = WENOZ5 erf.dryscal_vert_adv_type = WENOZ5 erf.moistscal_horiz_adv_type = WENOZ5 @@ -59,16 +61,25 @@ erf.moistscal_vert_adv_type = WENOZ5 erf.moisture_model = "SAM" erf.molec_diff_type = "None" + erf.les_type = "Smagorinsky" -erf.Cs = 0.1 +erf.Cs = 0.17 + +#erf.les_type = "Deardorff" +#erf.Ck = 0.1 +#erf.sigma_k = 1.0 +#erf.Ce = 0.1 + +erf.Pr_t = 0.333333 +erf.Sc_t = 0.333333 erf.init_type = "input_sounding" erf.init_sounding_ideal = true -erf.add_custom_rhotheta_forcing = true -erf.add_custom_moisture_forcing = true -erf.add_custom_geostrophic_profile = true -erf.add_custom_w_subsidence = true +erf.add_custom_rhotheta_forcing = true +erf.add_custom_moisture_forcing = true +erf.add_custom_geostrophic_profile = true +erf.add_custom_w_subsidence = true erf.custom_forcing_uses_primitive_vars = false # Higher values of perturbations lead to instability diff --git a/Exec/DevTests/Bomex/prob.H b/Exec/DevTests/Bomex/prob.H index ec14f6f0d..44ee89000 100644 --- a/Exec/DevTests/Bomex/prob.H +++ b/Exec/DevTests/Bomex/prob.H @@ -11,7 +11,7 @@ struct ProbParm : ProbParmDefaults { amrex::Real rho_0 = 0.0; amrex::Real T_0 = 0.0; amrex::Real A_0 = 1.0; - amrex::Real QKE_0 = 0.1; + amrex::Real KE_0 = 0.1; amrex::Real U_0 = 0.0; amrex::Real V_0 = 0.0; diff --git a/Exec/DevTests/Bomex/prob.cpp b/Exec/DevTests/Bomex/prob.cpp index 307979c05..80e46dd34 100644 --- a/Exec/DevTests/Bomex/prob.cpp +++ b/Exec/DevTests/Bomex/prob.cpp @@ -17,7 +17,7 @@ Problem::Problem (const Real* problo, const Real* probhi) pp.query("rho_0", parms.rho_0); pp.query("T_0", parms.T_0); pp.query("A_0", parms.A_0); - pp.query("QKE_0", parms.QKE_0); + pp.query("QKE_0", parms.KE_0); pp.query("U_0", parms.U_0); pp.query("V_0", parms.V_0); @@ -71,7 +71,7 @@ Problem::init_custom_pert ( Array4 const& x_vel_pert, Array4 const& y_vel_pert, Array4 const& z_vel_pert, - Array4 const& /*r_hse*/, + Array4 const& r_hse, Array4 const& /*p_hse*/, Array4 const& /*z_nd*/, Array4 const& /*z_cc*/, @@ -111,9 +111,10 @@ Problem::init_custom_pert ( // Set an initial value for QKE if (parms.custom_TKE) { - state_pert(i, j, k, RhoQKE_comp) = 1.0 - zc/3000.0; //*state_pert(i, j, k, Rho_comp); + state_pert(i, j, k, RhoKE_comp) = (1.0 - z/prob_hi[2]) * r_hse(i,j,k); + } else { + state_pert(i, j, k, RhoKE_comp) = parms.KE_0; } - else state_pert(i, j, k, RhoQKE_comp) = parms.QKE_0; if (use_moisture) { state_pert(i, j, k, RhoQ1_comp) = 0.0; diff --git a/Exec/EWP/inputs_1WT_lat_lon b/Exec/EWP/inputs_1WT_lat_lon new file mode 100644 index 000000000..fe77134fb --- /dev/null +++ b/Exec/EWP/inputs_1WT_lat_lon @@ -0,0 +1,106 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 3000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 100000.0 100000.0 1000 +amr.n_cell = 50 50 80 + +# WINDFARM PARAMETRIZATION PARAMETERS +erf.windfarm_type = "Fitch" +erf.windfarm_loc_type = "lat_lon" +erf.latitude_lo = 35.0 +erf.longitude_lo = -100.0 +erf.windfarm_loc_table = "windturbines_loc_lat_lon_1WT.txt" +erf.windfarm_spec_table = "windturbines_spec_1WT.tbl" + +#erf.grid_stretching_ratio = 1.025 +#erf.initial_dz = 16.0 + +geometry.is_periodic = 0 0 0 + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +#zlo.type = "MOST" +#erf.most.z0 = 0.1 +#erf.most.zref = 8.0 + +zlo.type = "SlipWall" +zhi.type = "SlipWall" +xlo.type = "Inflow" +xhi.type = "Outflow" +ylo.type = "Outflow" +yhi.type = "Outflow" + +xlo.velocity = 10. 0. 0. +xlo.density = 1.226 +xlo.theta = 300. + +#erf.sponge_strength = 0.1 +#erf.use_xlo_sponge_damping = true +#erf.xlo_sponge_end = 10000.0 +#erf.use_xhi_sponge_damping = true +#erf.xhi_sponge_start = 90000.0 + +#erf.sponge_density = 1.226 +#erf.sponge_x_velocity = 10.0 +#erf.sponge_y_velocity = 0.0 +#erf.sponge_z_velocity = 0.0 + + +# TIME STEP CONTROL +erf.use_native_mri = 1 +erf.fixed_dt = 3.0 # fixed time step depending on grid resolution +#erf.fixed_fast_dt = 0.0025 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 1000 # number of timesteps between checkpoints +#erf.restart = chk02000 + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta QKE num_turb vorticity + +# ADVECTION SCHEMES +erf.dycore_horiz_adv_type = "Centered_2nd" +erf.dycore_vert_adv_type = "Centered_2nd" +erf.dryscal_horiz_adv_type = "Centered_2nd" +erf.dryscal_vert_adv_type = "Centered_2nd" +erf.moistscal_horiz_adv_type = "Centered_2nd" +erf.moistscal_vert_adv_type = "Centered_2nd" + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "ConstantAlpha" +erf.les_type = "None" +erf.Cs = 1.5 +erf.dynamicViscosity = 10.0 + +erf.pbl_type = "None" + +erf.init_type = "uniform" + + +# PROBLEM PARAMETERS +prob.rho_0 = 1.226 +prob.A_0 = 1.0 + +prob.U_0 = 10.0 +prob.V_0 = 0.0 +prob.W_0 = 0.0 +prob.T_0 = 300.0 diff --git a/Exec/EWP/inputs_1WT_x_y b/Exec/EWP/inputs_1WT_x_y new file mode 100644 index 000000000..ee14f42cf --- /dev/null +++ b/Exec/EWP/inputs_1WT_x_y @@ -0,0 +1,104 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 3000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 100000.0 100000.0 1000 +amr.n_cell = 50 50 80 + +# WINDFARM PARAMETRIZATION PARAMETERS +erf.windfarm_type = "EWP" +erf.windfarm_loc_type = "x_y" +erf.windfarm_loc_table = "windturbines_loc_x_y_1WT.txt" +erf.windfarm_spec_table = "windturbines_spec_1WT.tbl" + +#erf.grid_stretching_ratio = 1.025 +#erf.initial_dz = 16.0 + +geometry.is_periodic = 0 0 0 + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +#zlo.type = "MOST" +#erf.most.z0 = 0.1 +#erf.most.zref = 8.0 + +zlo.type = "SlipWall" +zhi.type = "SlipWall" +xlo.type = "Inflow" +xhi.type = "Outflow" +ylo.type = "Outflow" +yhi.type = "Outflow" + +xlo.velocity = 10. 0. 0. +xlo.density = 1.226 +xlo.theta = 300. + +#erf.sponge_strength = 0.1 +#erf.use_xlo_sponge_damping = true +#erf.xlo_sponge_end = 10000.0 +#erf.use_xhi_sponge_damping = true +#erf.xhi_sponge_start = 90000.0 + +#erf.sponge_density = 1.226 +#erf.sponge_x_velocity = 10.0 +#erf.sponge_y_velocity = 0.0 +#erf.sponge_z_velocity = 0.0 + + +# TIME STEP CONTROL +erf.use_native_mri = 1 +erf.fixed_dt = 3.0 # fixed time step depending on grid resolution +#erf.fixed_fast_dt = 0.0025 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 1000 # number of timesteps between checkpoints +#erf.restart = chk02000 + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta QKE num_turb vorticity + +# ADVECTION SCHEMES +erf.dycore_horiz_adv_type = "Centered_2nd" +erf.dycore_vert_adv_type = "Centered_2nd" +erf.dryscal_horiz_adv_type = "Centered_2nd" +erf.dryscal_vert_adv_type = "Centered_2nd" +erf.moistscal_horiz_adv_type = "Centered_2nd" +erf.moistscal_vert_adv_type = "Centered_2nd" + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "ConstantAlpha" +erf.les_type = "None" +erf.Cs = 1.5 +erf.dynamicViscosity = 10.0 + +erf.pbl_type = "None" + +erf.init_type = "uniform" + + +# PROBLEM PARAMETERS +prob.rho_0 = 1.226 +prob.A_0 = 1.0 + +prob.U_0 = 10.0 +prob.V_0 = 0.0 +prob.W_0 = 0.0 +prob.T_0 = 300.0 diff --git a/Exec/EWP/inputs_WindFarm_lat_lon b/Exec/EWP/inputs_WindFarm_lat_lon new file mode 100644 index 000000000..e6ebae87b --- /dev/null +++ b/Exec/EWP/inputs_WindFarm_lat_lon @@ -0,0 +1,85 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 100000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 200150 202637 1000 +amr.n_cell = 50 50 40 + +# WINDFARM PARAMETRIZATION PARAMETERS +erf.windfarm_type = "Fitch" +erf.windfarm_loc_type = "lat_lon" +erf.latitude_lo = 35.0 +erf.longitude_lo = -100.0 +erf.windfarm_loc_table = "windturbines_loc_lat_lon_WindFarm.txt" +erf.windfarm_spec_table = "windturbines_spec_WindFarm.tbl" + + +#erf.grid_stretching_ratio = 1.025 +#erf.initial_dz = 16.0 + +geometry.is_periodic = 0 0 0 + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +#zlo.type = "MOST" +#erf.most.z0 = 0.1 +#erf.most.zref = 8.0 + +zlo.type = "SlipWall" +zhi.type = "SlipWall" +xlo.type = "Outflow" +xhi.type = "Outflow" +ylo.type = "Outflow" +yhi.type = "Outflow" + +# TIME STEP CONTROL +erf.use_native_mri = 1 +erf.fixed_dt = 0.25 # fixed time step depending on grid resolution +#erf.fixed_fast_dt = 0.0025 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 10000 # number of timesteps between checkpoints +#erf.restart = chk02000 + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 1000 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta QKE num_turb vorticity + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "ConstantAlpha" +erf.les_type = "None" +erf.Cs = 1.5 +erf.dynamicViscosity = 100.0 + +erf.pbl_type = "None" + +erf.init_type = "uniform" + +erf.windfarm_type = "EWP" + +# PROBLEM PARAMETERS +prob.rho_0 = 1.0 +prob.A_0 = 1.0 + +prob.U_0 = 10.0 +prob.V_0 = 10.0 +prob.W_0 = 0.0 +prob.T_0 = 300.0 + diff --git a/Exec/EWP/inputs_EWP_WindFarm b/Exec/EWP/inputs_WindFarm_x_y similarity index 89% rename from Exec/EWP/inputs_EWP_WindFarm rename to Exec/EWP/inputs_WindFarm_x_y index b7284ee28..7fbed1830 100644 --- a/Exec/EWP/inputs_EWP_WindFarm +++ b/Exec/EWP/inputs_WindFarm_x_y @@ -6,11 +6,16 @@ amrex.fpe_trap_invalid = 1 fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY -erf.latitude_lo = 35.0 -erf.longitude_lo = -100.0 geometry.prob_extent = 200150 202637 1000 amr.n_cell = 50 50 40 +# WINDFARM PARAMETRIZATION PARAMETERS +erf.windfarm_type = "Fitch" +erf.windfarm_loc_type = "x_y" +erf.windfarm_loc_table = "windturbines_loc_x_y_WindFarm.txt" +erf.windfarm_spec_table = "windturbines_spec_WindFarm.tbl" + + #erf.grid_stretching_ratio = 1.025 #erf.initial_dz = 16.0 diff --git a/Exec/EWP/windturbines.txt_LargeDomain b/Exec/EWP/windturbines.txt_LargeDomain deleted file mode 100644 index 586e94b81..000000000 --- a/Exec/EWP/windturbines.txt_LargeDomain +++ /dev/null @@ -1,30 +0,0 @@ -36.732 -97.128 1 -36.824 -97.367 1 -36.967 -96.754 1 -36.607 -96.846 1 -36.719 -97.249 1 -36.889 -97.031 1 -36.621 -96.932 1 -36.752 -96.788 1 -36.855 -97.421 1 -36.694 -96.656 1 -36.579 -97.223 1 -36.912 -97.379 1 -36.721 -96.815 1 -36.685 -97.437 1 -36.781 -97.049 1 -36.936 -97.203 1 -36.541 -96.791 1 -36.655 -97.284 1 -36.798 -96.695 1 -36.887 -96.956 1 -36.633 -96.716 1 -36.912 -96.813 1 -36.725 -96.579 1 -36.594 -97.184 1 -36.902 -97.354 1 -36.628 -96.871 1 -36.762 -97.428 1 -36.811 -97.285 1 -36.549 -97.492 1 -36.642 -96.673 1 diff --git a/Exec/EWP/windturbines_WindFarm.txt b/Exec/EWP/windturbines_WindFarm.txt deleted file mode 100644 index 4c68601f6..000000000 --- a/Exec/EWP/windturbines_WindFarm.txt +++ /dev/null @@ -1,97 +0,0 @@ -35.7828828829 -99.0168 1 -35.8219219219 -99.0168 1 -35.860960961 -99.0168 1 -35.9 -99.0168 1 -35.939039039 -99.0168 1 -35.9780780781 -99.0168 1 -36.0171171171 -99.0168 1 -35.7828828829 -98.9705333333 1 -35.8219219219 -98.9705333333 1 -35.860960961 -98.9705333333 1 -35.9 -98.9705333333 1 -35.939039039 -98.9705333333 1 -35.9780780781 -98.9705333333 1 -36.0171171171 -98.9705333333 1 -35.7828828829 -98.9242666667 1 -35.8219219219 -98.9242666667 1 -35.860960961 -98.9242666667 1 -35.9 -98.9242666667 1 -35.939039039 -98.9242666667 1 -35.9780780781 -98.9242666667 1 -36.0171171171 -98.9242666667 1 -35.7828828829 -98.878 1 -35.8219219219 -98.878 1 -35.860960961 -98.878 1 -35.9 -98.878 1 -35.939039039 -98.878 1 -35.9780780781 -98.878 1 -36.0171171171 -98.878 1 -35.7828828829 -98.8317333333 1 -35.8219219219 -98.8317333333 1 -35.860960961 -98.8317333333 1 -35.9 -98.8317333333 1 -35.939039039 -98.8317333333 1 -35.9780780781 -98.8317333333 1 -36.0171171171 -98.8317333333 1 -35.7828828829 -98.7854666667 1 -35.8219219219 -98.7854666667 1 -35.860960961 -98.7854666667 1 -35.9 -98.7854666667 1 -35.939039039 -98.7854666667 1 -35.9780780781 -98.7854666667 1 -36.0171171171 -98.7854666667 1 -35.7828828829 -98.7392 1 -35.8219219219 -98.7392 1 -35.860960961 -98.7392 1 -35.9 -98.7392 1 -35.939039039 -98.7392 1 -35.9780780781 -98.7392 1 -36.0171171171 -98.7392 1 -35.463963964 -98.2228 1 -35.487987988 -98.2228 1 -35.512012012 -98.2228 1 -35.536036036 -98.2228 1 -35.463963964 -98.1929333333 1 -35.487987988 -98.1929333333 1 -35.512012012 -98.1929333333 1 -35.536036036 -98.1929333333 1 -35.463963964 -98.1630666667 1 -35.487987988 -98.1630666667 1 -35.512012012 -98.1630666667 1 -35.536036036 -98.1630666667 1 -35.463963964 -98.1332 1 -35.487987988 -98.1332 1 -35.512012012 -98.1332 1 -35.536036036 -98.1332 1 -36.363963964 -99.4228 1 -36.387987988 -99.4228 1 -36.412012012 -99.4228 1 -36.436036036 -99.4228 1 -36.363963964 -99.3929333333 1 -36.387987988 -99.3929333333 1 -36.412012012 -99.3929333333 1 -36.436036036 -99.3929333333 1 -36.363963964 -99.3630666667 1 -36.387987988 -99.3630666667 1 -36.412012012 -99.3630666667 1 -36.436036036 -99.3630666667 1 -36.363963964 -99.3332 1 -36.387987988 -99.3332 1 -36.412012012 -99.3332 1 -36.436036036 -99.3332 1 -35.263963964 -99.5228 1 -35.287987988 -99.5228 1 -35.312012012 -99.5228 1 -35.336036036 -99.5228 1 -35.263963964 -99.4929333333 1 -35.287987988 -99.4929333333 1 -35.312012012 -99.4929333333 1 -35.336036036 -99.4929333333 1 -35.263963964 -99.4630666667 1 -35.287987988 -99.4630666667 1 -35.312012012 -99.4630666667 1 -35.336036036 -99.4630666667 1 -35.263963964 -99.4332 1 -35.287987988 -99.4332 1 -35.312012012 -99.4332 1 -35.336036036 -99.4332 1 diff --git a/Exec/EWP/windturbines_loc_lat_lon_1WT.txt b/Exec/EWP/windturbines_loc_lat_lon_1WT.txt new file mode 100644 index 000000000..b78b36903 --- /dev/null +++ b/Exec/EWP/windturbines_loc_lat_lon_1WT.txt @@ -0,0 +1,2 @@ +35.45 -99.5 1 + diff --git a/Exec/EWP/windturbines.txt b/Exec/EWP/windturbines_loc_lat_lon_WindFarm.txt similarity index 100% rename from Exec/EWP/windturbines.txt rename to Exec/EWP/windturbines_loc_lat_lon_WindFarm.txt diff --git a/Exec/EWP/windturbines_loc_x_y_1WT.txt b/Exec/EWP/windturbines_loc_x_y_1WT.txt new file mode 100644 index 000000000..e58d4be2e --- /dev/null +++ b/Exec/EWP/windturbines_loc_x_y_1WT.txt @@ -0,0 +1 @@ +45513.3107451298 49950.0000000004 diff --git a/Exec/EWP/windturbines_loc_x_y_WindFarm.txt b/Exec/EWP/windturbines_loc_x_y_WindFarm.txt new file mode 100644 index 000000000..19fe95746 --- /dev/null +++ b/Exec/EWP/windturbines_loc_x_y_WindFarm.txt @@ -0,0 +1,96 @@ +89264.99080053 91233.3333309002 +89259.1966417755 95566.6666710007 +89254.1277665419 99900.0000000001 +89249.7842982733 104233.333329 +89246.1663427532 108566.6666691 +89243.2739881117 112899.9999981 +93458.6633652711 86900.0000019001 +93450.4377452458 91233.3333309002 +93442.9032518779 95566.6666710007 +93436.0600522829 99900.0000000001 +93429.9082982192 104233.333329 +93424.4481261366 108566.6666691 +93419.6796571948 112899.9999981 +97646.3846285663 86900.0000019001 +97636.5111258549 91233.3333309002 +97627.2975408011 95566.6666710007 +97618.7440601867 99900.0000000001 +97610.850857392 104233.333329 +97603.6180924568 108566.6666691 +97597.045912112 112899.9999981 +101834.602975817 86900.0000019001 +101823.132838651 91233.3333309002 +101812.293876018 95566.6666710007 +101802.086289453 99900.0000000001 +101792.510268736 104233.333329 +101783.565991962 108566.6666691 +101775.253625592 112899.9999981 +106023.258628483 86900.0000019001 +106010.237048028 91233.3333309002 +105997.820076828 95566.6666710007 +105986.007927316 99900.0000000001 +105974.800801568 104233.333329 +105964.198891379 108566.6666691 +105954.202378331 112899.9999981 +110212.300853635 86900.0000019001 +110197.767881144 91233.3333309002 +110183.814885453 95566.6666710007 +110170.442086857 99900.0000000001 +110157.64969648 104233.333329 +110145.437916369 108566.6666691 +110133.806939566 112899.9999981 +114401.686336999 86900.0000019001 +114385.677633874 91233.3333309002 +114370.225998194 95566.6666710007 +114355.331655708 99900.0000000001 +114340.994824013 104233.333329 +114327.215712655 108566.6666691 +114313.994523223 112899.9999981 +161442.612044564 51500.0000039998 +161421.662101279 54166.666668 +161400.854802928 56833.3333320002 +161380.190204659 59499.9999960004 +164154.692973227 51500.0000039998 +164133.28852851 54166.666668 +164112.023937775 56833.3333320002 +164090.899255367 59499.9999960004 +166866.785694878 51500.0000039998 +166844.928396839 54166.666668 +166823.208245294 56833.3333320002 +166801.625293788 59499.9999960004 +169578.889421966 51500.0000039998 +169556.58083956 54166.666668 +169534.406775621 56833.3333320002 +169512.367282909 59499.9999960004 +52896.6844501924 151400.000004 +52915.9012759738 154066.666668 +52935.5802920911 156733.333332 +52955.7209832526 159399.999996 +55556.8294903758 151400.000004 +55574.333630956 154066.666668 +55592.2786959773 156733.333332 +55610.6642585664 159399.999996 +58220.3999614163 151400.000004 +58236.3087392296 154066.666668 +58252.6389095071 156733.333332 +58269.390117932 159399.999996 +60886.9461031228 151400.000004 +60901.3618319475 154066.666668 +60916.1809602466 156733.333332 +60931.4031936588 159399.999996 +43430.4219718844 29300.0000039997 +43430.6186514159 31966.666668 +43431.3883478816 34633.3333320001 +43432.7310308126 37299.9999960003 +46144.3886884569 29300.0000039997 +46143.7943466196 31966.666668 +46143.7390225352 34633.3333320001 +46144.2227181378 37299.9999960003 +48858.5817150023 29300.0000039997 +48857.2395489498 31966.666668 +48856.4061323542 34633.3333320001 +48856.0814912463 37299.9999960003 +51572.9651124459 29300.0000039997 +51570.9115091189 31966.666668 +51569.3395333094 34633.3333320001 +51568.2492290571 37299.9999960003 diff --git a/Exec/EWP/windturbines_spec_1WT.tbl b/Exec/EWP/windturbines_spec_1WT.tbl new file mode 100644 index 000000000..e1abbf4cf --- /dev/null +++ b/Exec/EWP/windturbines_spec_1WT.tbl @@ -0,0 +1,6 @@ +4 +119.0 178.0 0.130 2.0 +9 0.805 50.0 +10 0.805 50.0 +11 0.805 50.0 +12 0.805 50.0 diff --git a/Exec/EWP/wind-turbine-1.tbl b/Exec/EWP/windturbines_spec_WindFarm.tbl similarity index 100% rename from Exec/EWP/wind-turbine-1.tbl rename to Exec/EWP/windturbines_spec_WindFarm.tbl diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index 1e57be63f..855377953 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -39,6 +39,9 @@ public: amrex::Vector>>& lmask_lev, amrex::Vector> lsm_data, amrex::Vector> lsm_flux, + amrex::Vector>& Hwave, + amrex::Vector>& Lwave, + amrex::Vector>& eddyDiffs, amrex::Real start_bdy_time = 0.0, amrex::Real bdy_time_interval = 0.0) : m_exp_most(use_exp_most), @@ -117,6 +120,8 @@ public: } else if (rough_string == "modified_charnock") { rough_type = RoughCalcType::MODIFIED_CHARNOCK; pp.query("most.modified_charnock_depth",depth); + } else if (rough_string == "wave_coupled") { + rough_type = RoughCalcType::WAVE_COUPLED; } else { amrex::Abort("Undefined MOST roughness type!"); } @@ -156,6 +161,16 @@ public: } } + // Get pointers to wave data + m_Hwave_lev.resize(nlevs); + m_Lwave_lev.resize(nlevs); + m_eddyDiffs_lev.resize(nlevs); + for (int lev(0); lev @@ -235,7 +249,6 @@ public: amrex::MultiFab* xzmom_flux, amrex::MultiFab* zxmom_flux, amrex::MultiFab* yzmom_flux, amrex::MultiFab* zymom_flux, amrex::MultiFab* heat_flux, - amrex::MultiFab* eddyDiffs, amrex::MultiFab* z_phys, const amrex::Real& dz_no_terrain, const FluxCalc& flux_comp); @@ -306,7 +319,8 @@ public: enum struct RoughCalcType { CONSTANT = 0, ///< Constant z0 CHARNOCK, - MODIFIED_CHARNOCK + MODIFIED_CHARNOCK, + WAVE_COUPLED }; FluxCalcType flux_type{FluxCalcType::MOENG}; @@ -316,7 +330,7 @@ public: private: bool use_moisture; bool m_exp_most = false; - amrex::Real z0_const; + amrex::Real z0_const{0.1}; amrex::Real surf_temp; amrex::Real surf_heating_rate{0}; amrex::Real surf_temp_flux{0}; @@ -341,6 +355,9 @@ private: amrex::Vector> m_lmask_lev; amrex::Vector> m_lsm_data_lev; amrex::Vector> m_lsm_flux_lev; + amrex::Vector m_Hwave_lev; + amrex::Vector m_Lwave_lev; + amrex::Vector m_eddyDiffs_lev; }; #endif /* ABLMOST_H */ diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp index a9ebbcb7f..0f57c60c2 100644 --- a/Source/BoundaryConditions/ABLMost.cpp +++ b/Source/BoundaryConditions/ABLMost.cpp @@ -36,9 +36,12 @@ ABLMost::update_fluxes (const int& lev, } else if (rough_type == RoughCalcType::CHARNOCK) { surface_flux_charnock most_flux(m_ma.get_zref(), surf_temp_flux, cnk_a); compute_fluxes(lev, max_iters, most_flux); - } else { + } else if (rough_type == RoughCalcType::MODIFIED_CHARNOCK) { surface_flux_mod_charnock most_flux(m_ma.get_zref(), surf_temp_flux, depth); compute_fluxes(lev, max_iters, most_flux); + } else { + surface_flux_wave_coupled most_flux(m_ma.get_zref(), surf_temp_flux); + compute_fluxes(lev, max_iters, most_flux); } } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) { update_surf_temp(time); @@ -48,9 +51,12 @@ ABLMost::update_fluxes (const int& lev, } else if (rough_type == RoughCalcType::CHARNOCK) { surface_temp_charnock most_flux(m_ma.get_zref(), surf_temp_flux, cnk_a); compute_fluxes(lev, max_iters, most_flux); - } else { + } else if (rough_type == RoughCalcType::MODIFIED_CHARNOCK) { surface_temp_mod_charnock most_flux(m_ma.get_zref(), surf_temp_flux, depth); compute_fluxes(lev, max_iters, most_flux); + } else { + surface_temp_wave_coupled most_flux(m_ma.get_zref(), surf_temp_flux); + compute_fluxes(lev, max_iters, most_flux); } } else { if (rough_type == RoughCalcType::CONSTANT) { @@ -59,9 +65,12 @@ ABLMost::update_fluxes (const int& lev, } else if (rough_type == RoughCalcType::CHARNOCK) { adiabatic_charnock most_flux(m_ma.get_zref(), surf_temp_flux, cnk_a); compute_fluxes(lev, max_iters, most_flux); - } else { + } else if (rough_type == RoughCalcType::MODIFIED_CHARNOCK) { adiabatic_mod_charnock most_flux(m_ma.get_zref(), surf_temp_flux, depth); compute_fluxes(lev, max_iters, most_flux); + } else { + adiabatic_wave_coupled most_flux(m_ma.get_zref(), surf_temp_flux); + compute_fluxes(lev, max_iters, most_flux); } } // theta flux } else if (flux_type == FluxCalcType::CUSTOM) { @@ -102,10 +111,16 @@ ABLMost::compute_fluxes (const int& lev, const auto umm_arr = umm_ptr->array(mfi); const auto z0_arr = z_0[lev].array(); + // Wave properties if they exist + const auto Hwave_arr = (m_Hwave_lev[lev]) ? m_Hwave_lev[lev]->array(mfi) : Array4 {}; + const auto Lwave_arr = (m_Lwave_lev[lev]) ? m_Lwave_lev[lev]->array(mfi) : Array4 {}; + const auto eta_arr = m_eddyDiffs_lev[lev]->array(mfi); + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { most_flux.iterate_flux(i, j, k, max_iters, z0_arr, umm_arr, tm_arr, - u_star_arr, t_star_arr, t_surf_arr, olen_arr); + u_star_arr, t_star_arr, t_surf_arr, olen_arr, + Hwave_arr, Lwave_arr, eta_arr); }); } } @@ -124,7 +139,6 @@ ABLMost::impose_most_bcs (const int& lev, MultiFab* xzmom_flux, MultiFab* zxmom_flux, MultiFab* yzmom_flux, MultiFab* zymom_flux, MultiFab* heat_flux, - MultiFab* eddyDiffs, MultiFab* z_phys) { const int klo = 0; @@ -134,7 +148,6 @@ ABLMost::impose_most_bcs (const int& lev, xzmom_flux, zxmom_flux, yzmom_flux, zymom_flux, heat_flux, - eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp); } else if (flux_type == FluxCalcType::DONELAN) { donelan_flux flux_comp(klo); @@ -142,7 +155,6 @@ ABLMost::impose_most_bcs (const int& lev, xzmom_flux, zxmom_flux, yzmom_flux, zymom_flux, heat_flux, - eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp); } else { custom_flux flux_comp(klo); @@ -150,7 +162,6 @@ ABLMost::impose_most_bcs (const int& lev, xzmom_flux, zxmom_flux, yzmom_flux, zymom_flux, heat_flux, - eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp); } } @@ -171,7 +182,6 @@ ABLMost::compute_most_bcs (const int& lev, MultiFab* xzmom_flux, MultiFab* zxmom_flux, MultiFab* yzmom_flux, MultiFab* zymom_flux, MultiFab* heat_flux, - MultiFab* eddyDiffs, MultiFab* z_phys, const Real& dz_no_terrain, const FluxCalc& flux_comp) @@ -201,7 +211,7 @@ ABLMost::compute_most_bcs (const int& lev, auto t32_arr = (zymom_flux && m_exp_most) ? zymom_flux->array(mfi) : Array4{}; auto hfx_arr = (m_exp_most) ? heat_flux->array(mfi) : Array4{}; - const auto eta_arr = eddyDiffs->array(mfi); + const auto eta_arr = m_eddyDiffs_lev[lev]->array(mfi); const auto zphys_arr = (z_phys) ? z_phys->const_array(mfi) : Array4{}; diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index a63253297..dc630079b 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -386,7 +386,6 @@ ERF::FillIntermediatePatch (int lev, Real time, Tau13_lev[lev].get(), Tau31_lev[lev].get(), Tau23_lev[lev].get(), Tau32_lev[lev].get(), SFS_hfx3_lev[lev].get(), - eddyDiffs_lev[lev].get(), z_phys_nd[lev].get()); } diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index 2e4a45649..7f4931718 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -88,7 +88,10 @@ struct adiabatic const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, const amrex::Array4& /*t_surf_arr*/, - const amrex::Array4& olen_arr) const + const amrex::Array4& olen_arr, + const amrex::Array4& /*Hwave_arr*/, + const amrex::Array4& /*Lwave_arr*/, + const amrex::Array4& /*eta_arr*/) const { u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k)); t_star_arr(i,j,k) = 0.0; @@ -128,7 +131,10 @@ struct adiabatic_charnock const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, const amrex::Array4& /*t_surf_arr*/, - const amrex::Array4& olen_arr) const + const amrex::Array4& olen_arr, + const amrex::Array4& /*Hwave_arr*/, + const amrex::Array4& /*Lwave_arr*/, + const amrex::Array4& /*eta_arr*/) const { int iter = 0; amrex::Real ustar = 0.0; @@ -181,7 +187,10 @@ struct adiabatic_mod_charnock const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, const amrex::Array4& /*t_surf_arr*/, - const amrex::Array4& olen_arr) const + const amrex::Array4& olen_arr, + const amrex::Array4& /*Hwave_arr*/, + const amrex::Array4& /*Lwave_arr*/, + const amrex::Array4& /*eta_arr*/) const { int iter = 0; amrex::Real ustar = 0.0; @@ -206,6 +215,65 @@ private: }; +/** + * Adiabatic with wave-coupled roughness + */ +struct adiabatic_wave_coupled +{ + adiabatic_wave_coupled (amrex::Real zref, + amrex::Real flux) + { + mdata.zref = zref; + mdata.surf_temp_flux = flux; + } + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + void + iterate_flux (const int& i, + const int& j, + const int& k, + const int& max_iters, + const amrex::Array4& z0_arr, + const amrex::Array4& umm_arr, + const amrex::Array4& /*tm_arr*/, + const amrex::Array4& u_star_arr, + const amrex::Array4& t_star_arr, + const amrex::Array4& /*t_surf_arr*/, + const amrex::Array4& olen_arr, + const amrex::Array4& Hwave_arr, + const amrex::Array4& Lwave_arr, + const amrex::Array4& eta_arr) const + { + int iter = 0; + amrex::Real ustar = 0.0; + amrex::Real z0 = 0.0; + int ie, je; + ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; + je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k)); + do { + ustar = u_star_arr(i,j,k); + z0 = 1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/Lwave_arr(i,j,k), 4.5 ) + + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar; + u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0); + ++iter; + } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters); + + t_star_arr(i,j,k) = 0.0; + olen_arr(i,j,k) = 1.0e16; + z0_arr(i,j,k) = z0; + } + +private: + most_data mdata; + similarity_funs sfuns; + const amrex::Real tol = 1.0e-5; +}; + + /** * Surface flux with constant roughness */ @@ -231,7 +299,10 @@ struct surface_flux const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, const amrex::Array4& t_surf_arr, - const amrex::Array4& olen_arr) const + const amrex::Array4& olen_arr, + const amrex::Array4& /*Hwave_arr*/, + const amrex::Array4& /*Lwave_arr*/, + const amrex::Array4& /*eta_arr*/) const { int iter = 0; amrex::Real ustar = 0.0; @@ -291,7 +362,10 @@ struct surface_flux_charnock const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, const amrex::Array4& t_surf_arr, - const amrex::Array4& olen_arr) const + const amrex::Array4& olen_arr, + const amrex::Array4& /*Hwave_arr*/, + const amrex::Array4& /*Lwave_arr*/, + const amrex::Array4& /*eta_arr*/) const { int iter = 0; amrex::Real ustar = 0.0; @@ -355,7 +429,10 @@ struct surface_flux_mod_charnock const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, const amrex::Array4& t_surf_arr, - const amrex::Array4& olen_arr) const + const amrex::Array4& olen_arr, + const amrex::Array4& /*Hwave_arr*/, + const amrex::Array4& /*Lwave_arr*/, + const amrex::Array4& /*eta_arr*/) const { int iter = 0; amrex::Real ustar = 0.0; @@ -391,6 +468,76 @@ private: }; +/** + * Surface flux with wave-coupled roughness + */ +struct surface_flux_wave_coupled +{ + surface_flux_wave_coupled (amrex::Real zref, + amrex::Real flux) + { + mdata.zref = zref; + mdata.surf_temp_flux = flux; + } + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + void + iterate_flux (const int& i, + const int& j, + const int& k, + const int& max_iters, + const amrex::Array4& z0_arr, + const amrex::Array4& umm_arr, + const amrex::Array4& tm_arr, + const amrex::Array4& u_star_arr, + const amrex::Array4& t_star_arr, + const amrex::Array4& t_surf_arr, + const amrex::Array4& olen_arr, + const amrex::Array4& Hwave_arr, + const amrex::Array4& Lwave_arr, + const amrex::Array4& eta_arr) const + { + int iter = 0; + amrex::Real ustar = 0.0; + amrex::Real z0 = 0.0; + amrex::Real zeta = 0.0; + amrex::Real psi_m = 0.0; + amrex::Real psi_h = 0.0; + amrex::Real Olen = 0.0; + int ie, je; + ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; + je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k)); + do { + ustar = u_star_arr(i,j,k); + z0 = 1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/Lwave_arr(i,j,k), 4.5 ) + + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar; + Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / + (mdata.kappa * mdata.gravity * mdata.surf_temp_flux); + zeta = mdata.zref / Olen; + psi_m = sfuns.calc_psi_m(zeta); + psi_h = sfuns.calc_psi_h(zeta); + u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / (std::log(mdata.zref / z0) - psi_m); + ++iter; + } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters); + + t_surf_arr(i,j,k) = mdata.surf_temp_flux * (std::log(mdata.zref / z0) - psi_h) / + (u_star_arr(i,j,k) * mdata.kappa) + tm_arr(i,j,k); + t_star_arr(i,j,k) = -mdata.surf_temp_flux / u_star_arr(i,j,k); + olen_arr(i,j,k) = Olen; + z0_arr(i,j,k) = z0; + } + +private: + most_data mdata; + similarity_funs sfuns; + const amrex::Real tol = 1.0e-5; +}; + + /** * Surface temperature with constant roughness */ @@ -416,7 +563,10 @@ struct surface_temp const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, const amrex::Array4& t_surf_arr, - const amrex::Array4& olen_arr) const + const amrex::Array4& olen_arr, + const amrex::Array4& /*Hwave_arr*/, + const amrex::Array4& /*Lwave_arr*/, + const amrex::Array4& /*eta_arr*/) const { int iter = 0; amrex::Real ustar = 0.0; @@ -478,7 +628,10 @@ struct surface_temp_charnock const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, const amrex::Array4& t_surf_arr, - const amrex::Array4& olen_arr) const + const amrex::Array4& olen_arr, + const amrex::Array4& /*Hwave_arr*/, + const amrex::Array4& /*Lwave_arr*/, + const amrex::Array4& /*eta_arr*/) const { int iter = 0; amrex::Real ustar = 0.0; @@ -544,7 +697,10 @@ struct surface_temp_mod_charnock const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, const amrex::Array4& t_surf_arr, - const amrex::Array4& olen_arr) const + const amrex::Array4& olen_arr, + const amrex::Array4& /*Hwave_arr*/, + const amrex::Array4& /*Lwave_arr*/, + const amrex::Array4& /*eta_arr*/) const { int iter = 0; amrex::Real ustar = 0.0; @@ -582,6 +738,78 @@ private: }; +/** + * Surface temperature with wave-coupled roughness + */ +struct surface_temp_wave_coupled +{ + surface_temp_wave_coupled (amrex::Real zref, + amrex::Real flux) + { + mdata.zref = zref; + mdata.surf_temp_flux = flux; + } + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + void + iterate_flux (const int& i, + const int& j, + const int& k, + const int& max_iters, + const amrex::Array4& z0_arr, + const amrex::Array4& umm_arr, + const amrex::Array4& tm_arr, + const amrex::Array4& u_star_arr, + const amrex::Array4& t_star_arr, + const amrex::Array4& t_surf_arr, + const amrex::Array4& olen_arr, + const amrex::Array4& Hwave_arr, + const amrex::Array4& Lwave_arr, + const amrex::Array4& eta_arr) const + { + int iter = 0; + amrex::Real ustar = 0.0; + amrex::Real z0 = 0.0; + amrex::Real tflux = 0.0; + amrex::Real zeta = 0.0; + amrex::Real psi_m = 0.0; + amrex::Real psi_h = 0.0; + amrex::Real Olen = 0.0; + int ie, je; + ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; + je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k)); + do { + ustar = u_star_arr(i,j,k); + z0 = 1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/Lwave_arr(i,j,k), 4.5 ) + + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar; + tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa / + (std::log(mdata.zref / z0) - psi_h); + Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / + (mdata.kappa * mdata.gravity * tflux); + zeta = mdata.zref / Olen; + psi_m = sfuns.calc_psi_m(zeta); + psi_h = sfuns.calc_psi_h(zeta); + u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / (std::log(mdata.zref / z0) - psi_m); + ++iter; + } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters); + + t_star_arr(i,j,k) = mdata.kappa * (tm_arr(i,j,k) - t_surf_arr(i,j,k)) / + (std::log(mdata.zref / z0) - psi_h); + olen_arr(i,j,k) = Olen; + z0_arr(i,j,k) = z0; + } + +private: + most_data mdata; + similarity_funs sfuns; + const amrex::Real tol = 1.0e-5; +}; + + /** * Moeng flux formulation */ diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index a3351c756..f3ebcb0b4 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -40,6 +40,10 @@ enum struct WindFarmType { Fitch, EWP, None }; +enum struct WindFarmLocType{ + lat_lon, x_y, None +}; + enum struct LandSurfaceType { SLM, MM5, None }; @@ -306,8 +310,26 @@ struct SolverChoice { windfarm_type = WindFarmType::EWP; } else if (windfarm_type_string != "None") { - amrex::Abort("Dont know this windfarm_type. Currently only Fitch and EWP models are supported."); + amrex::Abort("Are you using windfarms? Dont know this windfarm_type. windfarm_type" + " has to be Fitch or EWP or None."); + } + + static std::string windfarm_loc_type_string = "None"; + windfarm_loc_type = WindFarmLocType::None; + pp.query("windfarm_loc_type", windfarm_loc_type_string); + if (windfarm_loc_type_string == "lat_lon") { + windfarm_loc_type = WindFarmLocType::lat_lon; } + else if (windfarm_loc_type_string == "x_y") { + windfarm_loc_type = WindFarmLocType::x_y; + } + else if (windfarm_loc_type_string != "None") { + amrex::Abort("Are you using windfarms? Dont know this windfarm_loc_type." + " windfarm_loc_type has to be specified as lat_lon or x_y."); + } + + pp.query("windfarm_loc_table", windfarm_loc_table); + pp.query("windfarm_spec_table", windfarm_spec_table); // Test if time averaged data is to be output pp.query("time_avg_vel",time_avg_vel); @@ -473,6 +495,7 @@ struct SolverChoice { TerrainType terrain_type; MoistureType moisture_type; WindFarmType windfarm_type; + WindFarmLocType windfarm_loc_type; LandSurfaceType lsm_type; ABLDriverType abl_driver_type; @@ -485,6 +508,7 @@ struct SolverChoice { bool do_precip {true}; bool use_moist_background {false}; - amrex::Real latitude_lo, longitude_lo; + amrex::Real latitude_lo=-1e10, longitude_lo=-1e10; + std::string windfarm_loc_table, windfarm_spec_table; }; #endif diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index db7742162..cf681c12a 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -110,101 +110,101 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, //*********************************************************************************** else if (turbChoice.les_type == LESType::Deardorff) { - const Real l_C_k = turbChoice.Ck; - const Real l_C_e = turbChoice.Ce; - const Real l_C_e_wall = turbChoice.Ce_wall; - const Real Ce_lcoeff = amrex::max(0.0, l_C_e - 1.9*l_C_k); - const Real l_abs_g = const_grav; - const Real l_inv_theta0 = 1.0 / turbChoice.theta_ref; + const Real l_C_k = turbChoice.Ck; + const Real l_C_e = turbChoice.Ce; + const Real l_C_e_wall = turbChoice.Ce_wall; + const Real Ce_lcoeff = amrex::max(0.0, l_C_e - 1.9*l_C_k); + const Real l_abs_g = const_grav; + const Real l_inv_theta0 = 1.0 / turbChoice.theta_ref; #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box bxcc = mfi.tilebox(); + for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box bxcc = mfi.tilebox(); - const Array4& mu_turb = eddyViscosity.array(mfi); - const Array4& hfx_x = Hfx1.array(mfi); - const Array4& hfx_y = Hfx2.array(mfi); - const Array4& hfx_z = Hfx3.array(mfi); - const Array4& diss = Diss.array(mfi); + const Array4& mu_turb = eddyViscosity.array(mfi); + const Array4& hfx_x = Hfx1.array(mfi); + const Array4& hfx_y = Hfx2.array(mfi); + const Array4& hfx_z = Hfx3.array(mfi); + const Array4& diss = Diss.array(mfi); - const Array4 &cell_data = cons_in.array(mfi); + const Array4 &cell_data = cons_in.array(mfi); - Array4 mf_u = mapfac_u.array(mfi); - Array4 mf_v = mapfac_v.array(mfi); + Array4 mf_u = mapfac_u.array(mfi); + Array4 mf_v = mapfac_v.array(mfi); - Array4 z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) : Array4{}; + Array4 z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) : Array4{}; - ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real dxInv = cellSizeInv[0]; - Real dyInv = cellSizeInv[1]; - Real dzInv = cellSizeInv[2]; - if (use_terrain) { - // the terrain grid is only deformed in z for now - dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr); - } - Real cellVolMsf = 1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv); - Real DeltaMsf = std::pow(cellVolMsf,1.0/3.0); - - // Calculate stratification-dependent mixing length (Deardorff 1980) - Real eps = std::numeric_limits::epsilon(); - Real dtheta_dz; - if (use_most && k==klo) { - if (exp_most) { - dtheta_dz = ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) - - cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) )*dzInv; - } else { - dtheta_dz = 0.5 * (-3 * cell_data(i,j,k ,RhoTheta_comp) - / cell_data(i,j,k ,Rho_comp) - + 4 * cell_data(i,j,k+1,RhoTheta_comp) - / cell_data(i,j,k+1,Rho_comp) - - cell_data(i,j,k+2,RhoTheta_comp) - / cell_data(i,j,k+2,Rho_comp) ) * dzInv; - } - } else { - dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) - - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv; - } - Real E = cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp); - Real strat = l_abs_g * dtheta_dz * l_inv_theta0; // stratification - Real length; - if (strat <= eps) { - length = DeltaMsf; - } else { - length = 0.76 * std::sqrt(E / strat); - // mixing length should be _reduced_ for stable stratification - length = amrex::min(length, DeltaMsf); - // following WRF, make sure the mixing length isn't too small - length = amrex::max(length, 0.001 * DeltaMsf); - } + ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real dxInv = cellSizeInv[0]; + Real dyInv = cellSizeInv[1]; + Real dzInv = cellSizeInv[2]; + if (use_terrain) { + // the terrain grid is only deformed in z for now + dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr); + } + Real cellVolMsf = 1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv); + Real DeltaMsf = std::pow(cellVolMsf,1.0/3.0); + + // Calculate stratification-dependent mixing length (Deardorff 1980) + Real eps = std::numeric_limits::epsilon(); + Real dtheta_dz; + if (use_most && k==klo) { + if (exp_most) { + dtheta_dz = ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) + - cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) )*dzInv; + } else { + dtheta_dz = 0.5 * (-3 * cell_data(i,j,k ,RhoTheta_comp) + / cell_data(i,j,k ,Rho_comp) + + 4 * cell_data(i,j,k+1,RhoTheta_comp) + / cell_data(i,j,k+1,Rho_comp) + - cell_data(i,j,k+2,RhoTheta_comp) + / cell_data(i,j,k+2,Rho_comp) ) * dzInv; + } + } else { + dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) + - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv; + } + Real E = cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp); + Real strat = l_abs_g * dtheta_dz * l_inv_theta0; // stratification + Real length; + if (strat <= eps) { + length = DeltaMsf; + } else { + length = 0.76 * std::sqrt(E / strat); + // mixing length should be _reduced_ for stable stratification + length = amrex::min(length, DeltaMsf); + // following WRF, make sure the mixing length isn't too small + length = amrex::max(length, 0.001 * DeltaMsf); + } - // Calculate eddy diffusivities - // K = rho * C_k * l * KE^(1/2) - mu_turb(i,j,k,EddyDiff::Mom_h) = cell_data(i,j,k,Rho_comp) * l_C_k * length * std::sqrt(E); - mu_turb(i,j,k,EddyDiff::Mom_v) = mu_turb(i,j,k,EddyDiff::Mom_h); - // KH = (1 + 2*l/delta) * mu_turb - mu_turb(i,j,k,EddyDiff::Theta_v) = (1.+2.*length/DeltaMsf) * mu_turb(i,j,k,EddyDiff::Mom_v); - - // Calculate SFS quantities - // - dissipation - Real Ce; - if ((l_C_e_wall > 0) && (k==0)) { - Ce = l_C_e_wall; - } else { - Ce = 1.9*l_C_k + Ce_lcoeff*length / DeltaMsf; - } - diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,1.5) / length; - // - heat flux - // (Note: If using ERF_EXPLICIT_MOST_STRESS, the value at k=0 will - // be overwritten when BCs are applied) - hfx_x(i,j,k) = 0.0; - hfx_y(i,j,k) = 0.0; - hfx_z(i,j,k) = -mu_turb(i,j,k,EddyDiff::Theta_v) * dtheta_dz; // (rho*w)' theta' [kg m^-2 s^-1 K] - }); - } + // Calculate eddy diffusivities + // K = rho * C_k * l * KE^(1/2) + mu_turb(i,j,k,EddyDiff::Mom_h) = cell_data(i,j,k,Rho_comp) * l_C_k * length * std::sqrt(E); + mu_turb(i,j,k,EddyDiff::Mom_v) = mu_turb(i,j,k,EddyDiff::Mom_h); + // KH = (1 + 2*l/delta) * mu_turb + mu_turb(i,j,k,EddyDiff::Theta_v) = (1.+2.*length/DeltaMsf) * mu_turb(i,j,k,EddyDiff::Mom_v); + + // Calculate SFS quantities + // - dissipation + Real Ce; + if ((l_C_e_wall > 0) && (k==0)) { + Ce = l_C_e_wall; + } else { + Ce = 1.9*l_C_k + Ce_lcoeff*length / DeltaMsf; + } + diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,1.5) / length; + // - heat flux + // (Note: If using ERF_EXPLICIT_MOST_STRESS, the value at k=0 will + // be overwritten when BCs are applied) + hfx_x(i,j,k) = 0.0; + hfx_y(i,j,k) = 0.0; + hfx_z(i,j,k) = -mu_turb(i,j,k,EddyDiff::Theta_v) * dtheta_dz; // (rho*w)' theta' [kg m^-2 s^-1 K] + }); + } } // Extrapolate Kturb in x/y, fill remaining elements (relevent to lev==0) @@ -271,6 +271,38 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, }); } + // Copy Theta_v component into lateral ghost cells if using Deardorff (populated above) + if (use_KE) { + if (i_lo == domain.smallEnd(0)) { + ParallelFor(planex, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int lj = amrex::min(amrex::max(j, domain.smallEnd(1)), domain.bigEnd(1)); + mu_turb(i_lo-i, j, k, EddyDiff::Theta_v) = mu_turb(i_lo, lj, k, EddyDiff::Theta_v); + }); + } + if (i_hi == domain.bigEnd(0)) { + ParallelFor(planex, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int lj = amrex::min(amrex::max(j, domain.smallEnd(1)), domain.bigEnd(1)); + mu_turb(i_hi+i, j, k, EddyDiff::Theta_v) = mu_turb(i_hi, lj, k, EddyDiff::Theta_v); + }); + } + if (j_lo == domain.smallEnd(1)) { + ParallelFor(planey, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int li = amrex::min(amrex::max(i, domain.smallEnd(0)), domain.bigEnd(0)); + mu_turb(i, j_lo-j, k, EddyDiff::Theta_v) = mu_turb(li, j_lo, k, EddyDiff::Theta_v); + }); + } + if (j_hi == domain.bigEnd(1)) { + ParallelFor(planey, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int li = amrex::min(amrex::max(i, domain.smallEnd(0)), domain.bigEnd(0)); + mu_turb(i, j_hi+j, k, EddyDiff::Theta_v) = mu_turb(li, j_hi, k, EddyDiff::Theta_v); + }); + } + } + // refactor the code to eliminate the need for ifdef's for (auto n = 1; n < (EddyDiff::NumDiffs-1)/2; ++n) { int offset = (EddyDiff::NumDiffs-1)/2; @@ -281,10 +313,10 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, if(use_QKE) { ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n; - int indx_v = indx + offset; - mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; - mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); + int indx = n; + int indx_v = indx + offset; + mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; + mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); }); } break; @@ -292,20 +324,21 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, if (use_KE) { ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n; - int indx_v = indx + offset; - mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; - mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); + int indx = n; + int indx_v = indx + offset; + mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; + mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); }); } break; default: ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n; - int indx_v = indx + offset; - mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; - mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); + int indx = n; + int indx_v = indx + offset; + mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; + // NOTE: Theta_v has already been set for Deardorff + if (!(indx_v == EddyDiff::Theta_v && use_KE)) mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); }); break; } @@ -342,12 +375,12 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, if(use_QKE) { ParallelFor(planez, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n; - int indx_v = indx + offset; - mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); - mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); - mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); - mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); + int indx = n; + int indx_v = indx + offset; + mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); + mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); + mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); + mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); }); } break; @@ -355,24 +388,24 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, if (use_KE) { ParallelFor(planez, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n; - int indx_v = indx + offset; - mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); - mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); - mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); - mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); + int indx = n; + int indx_v = indx + offset; + mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); + mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); + mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); + mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); }); } break; default: ParallelFor(planez, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n ; - int indx_v = indx + offset; - mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); - mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); - mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); - mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); + int indx = n; + int indx_v = indx + offset; + mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); + mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); + mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); + mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); }); break; } diff --git a/Source/ERF.H b/Source/ERF.H index d3712f66d..e7122ed5a 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -366,6 +366,14 @@ public: #ifdef ERF_USE_WINDFARM void init_windfarm(int lev); + void init_windfarm_lat_lon(const int lev, + const std::string windfarm_loc_table, + amrex::Vector& xloc, + amrex::Vector& yloc); + void init_windfarm_x_y(const int lev, + const std::string windfarm_loc_table, + amrex::Vector& xloc, + amrex::Vector& yloc); #endif #ifdef ERF_USE_EB @@ -630,12 +638,11 @@ private: std::unique_ptr micro; amrex::Vector> qmoist; // (lev,ncomp) This has up to 8 components: qt, qv, qc, qi, qp, qr, qs, qg + // Variables for wind farm parametrization models + amrex::Vector Nturb; amrex::Vector vars_windfarm; // Fitch: Vabs, Vabsdt, dudt, dvdt, dTKEdt // EWP: dudt, dvdt, dTKEdt - amrex::Real hub_height, rotor_dia, thrust_coeff_standing, nominal_power; - amrex::Vector wind_speed, thrust_coeff, power; - LandSurface lsm; amrex::Vector> lsm_data; // (lev,ncomp) Components: theta, q1, q2 @@ -705,6 +712,10 @@ private: amrex::Vector base_state; amrex::Vector base_state_new; + // Wave coupling data + amrex::Vector> Hwave; + amrex::Vector> Lwave; + // array of flux registers amrex::Vector advflux_reg; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 68c267430..8a4b04d92 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -268,6 +268,15 @@ ERF::ERF () base_state.resize(nlevs_max); base_state_new.resize(nlevs_max); + // Wave coupling data + Hwave.resize(nlevs_max); + Lwave.resize(nlevs_max); + for (int lev = 0; lev < max_level; ++lev) + { + Hwave[lev] = nullptr; + Lwave[lev] = nullptr; + } + // Theta prim for MOST Theta_prim.resize(nlevs_max); @@ -910,7 +919,7 @@ ERF::InitData () } m_most = std::make_unique(geom, use_exp_most, vars_old, Theta_prim, Qv_prim, z_phys_nd, - sst_lev, lmask_lev, lsm_data, lsm_flux + sst_lev, lmask_lev, lsm_data, lsm_flux, Hwave, Lwave, eddyDiffs_lev #ifdef ERF_USE_NETCDF ,start_bdy_time, bdy_time_interval #endif @@ -1904,6 +1913,15 @@ ERF::ERF (const RealBox& rb, int max_level_in, base_state.resize(nlevs_max); base_state_new.resize(nlevs_max); + // Wave coupling data + Hwave.resize(nlevs_max); + Lwave.resize(nlevs_max); + for (int lev = 0; lev < max_level; ++lev) + { + Hwave[lev] = nullptr; + Lwave[lev] = nullptr; + } + // Theta prim for MOST Theta_prim.resize(nlevs_max); diff --git a/Source/Initialization/ERF_init_windfarm.cpp b/Source/Initialization/ERF_init_windfarm.cpp index 7963ac01d..fdaf02765 100644 --- a/Source/Initialization/ERF_init_windfarm.cpp +++ b/Source/Initialization/ERF_init_windfarm.cpp @@ -3,6 +3,7 @@ */ #include +#include using namespace amrex; @@ -15,39 +16,17 @@ using namespace amrex; void ERF::init_windfarm (int lev) { - // Read turbine locations from windturbines.txt - std::ifstream file("windturbines.txt"); - if (!file.is_open()) { - amrex::Error("Wind turbines location file windturbines.txt not found"); + Vector xloc, yloc; + if(solverChoice.windfarm_loc_type == WindFarmLocType::lat_lon) { + init_windfarm_lat_lon(lev, solverChoice.windfarm_loc_table, xloc, yloc); } - // Vector of vectors to store the matrix - std::vector lat, lon, xloc, yloc; - Real value1, value2, value3; - while (file >> value1 >> value2 >> value3) { - lat.push_back(value1); - lon.push_back(value2); + else if(solverChoice.windfarm_loc_type == WindFarmLocType::x_y) { + init_windfarm_x_y(lev, solverChoice.windfarm_loc_table, xloc, yloc); } - file.close(); - - Real rad_earth = 6371.0e3; // Radius of the earth - Real lat_lo = solverChoice.latitude_lo*M_PI/180.0; - Real lon_lo = solverChoice.longitude_lo*M_PI/180.0; - - // (lat_lo, lon_lo) is mapped to (0,0) - - for(int it=0;it& xloc, + Vector& yloc) +{ - //The first line is the number of pairs entries for the power curve and thrust coefficient. - //The second line gives first the height in meters of the turbine hub, second, the diameter in - //meters of the rotor, third the standing thrust coefficient, and fourth the nominal power of - //the turbine in MW. - //The remaining lines contain the three values of: wind speed, thrust coefficient, and power production in kW. + if(solverChoice.latitude_lo == -1e10) { + amrex::Error("Are you using latitude-longitude for wind turbine locations? There needs to be" + " entry for erf.latitude_lo in the inputs for" + " the lower bottom corner of the domain"); + } - // Read turbine data from wind-turbine-1.tbl - std::ifstream file_turb_table("wind-turbine-1.tbl"); - if (!file_turb_table.is_open()) { - amrex::Error("Wind turbines location file wind-turbine-1.tbl not found"); + if(solverChoice.longitude_lo == -1e10) { + amrex::Error("Are you using latitude-longitude for wind turbine locations? There needs to be" + " entry erf.longitude_lo in the inputs for" + " the lower bottom corner of the domain"); } - int nlines; - file_turb_table >> nlines; - wind_speed.resize(nlines); - thrust_coeff.resize(nlines); - power.resize(nlines); + if(std::fabs(solverChoice.latitude_lo) > 90.0) { + amrex::Error("The value of erf.latitude_lo in the input file should be within -90 and 90"); + } - file_turb_table >> hub_height >> rotor_dia >> thrust_coeff_standing >> nominal_power; - if(rotor_dia/2.0 > hub_height) - { - amrex::Abort("The blade length is more than the hub height. Check the second line in wind-turbine-1.tbl. Aborting....."); + if(std::fabs(solverChoice.longitude_lo) > 180.0) { + amrex::Error("The value of erf.longitude_lo in the input file should be within -180 and 180"); } - for(int iline=0;iline> wind_speed[iline] >> thrust_coeff[iline] >> power[iline]; + // Read turbine locations from windturbines.txt + std::ifstream file(windfarm_loc_table); + if (!file.is_open()) { + amrex::Error("Wind turbines location table not found. Either the inputs is missing the" + " erf.windfarm_loc_table entry or the file specified in the entry " + windfarm_loc_table + " is missing."); + } + // Vector of vectors to store the matrix + Vector lat, lon; + Real value1, value2, value3; + + while (file >> value1 >> value2 >> value3) { + + if(std::fabs(value1) > 90.0) { + amrex::Error("The value of latitude for entry " + std::to_string(lat.size() + 1) + + " in " + windfarm_loc_table + " should be within -90 and 90"); + } + + if(std::fabs(value2) > 180.0) { + amrex::Error("The value of longitude for entry " + std::to_string(lat.size() + 1) + + " in " + windfarm_loc_table + " should be within -180 and 180"); + } + + if(std::fabs(solverChoice.longitude_lo) > 180.0) { + amrex::Error("The values of erf.longitude_lo should be within -180 and 180"); + } + lat.push_back(value1); + lon.push_back(value2); } file.close(); + + Real rad_earth = 6371.0e3; // Radius of the earth + + Real lat_lo = solverChoice.latitude_lo*M_PI/180.0; + Real lon_lo = solverChoice.longitude_lo*M_PI/180.0; + + // (lat_lo, lon_lo) is mapped to (0,0) + + for(int it=0;it& xloc, + Vector& yloc) +{ + // Read turbine locations from windturbines.txt + std::ifstream file(windfarm_loc_table); + if (!file.is_open()) { + amrex::Error("Wind turbines location table not found. Either the inputs is missing the" + " erf.windfarm_loc_table entry or the file specified in the entry " + windfarm_loc_table + " is missing."); + } + // Vector of vectors to store the matrix + Real value1, value2; + while (file >> value1 >> value2) { + xloc.push_back(value1); + yloc.push_back(value2); + } + file.close(); +} diff --git a/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp b/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp index fbe9f506e..d92f3ded9 100644 --- a/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp +++ b/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp @@ -1,10 +1,12 @@ +#include +#include #include #include + using namespace amrex; -Real R_ewp = 30.0; -Real z_c_ewp = 100.0; -Real C_T_ewp = 0.8, C_TKE_ewp = 0.0; +Real C_TKE_ewp = 0.0; +Real K_turb = 1.0; void ewp_advance (int lev, const Geometry& geom, @@ -60,6 +62,7 @@ void ewp_source_terms_cellcentered (const Geometry& geom, auto dx = geom.CellSizeArray(); auto ProbHiArr = geom.ProbHiArray(); auto ProbLoArr = geom.ProbLoArray(); + Real sigma_0 = 1.7*rotor_rad; // Domain valid box const amrex::Box& domain = geom.Domain(); @@ -102,13 +105,13 @@ void ewp_source_terms_cellcentered (const Geometry& geom, v_vel(i,j,k)*v_vel(i,j,k) + w_vel(i,j,kk)*w_vel(i,j,kk), 0.5); - Real L_wake = std::pow(dx[0]*dx[1],0.5); - Real K_turb = 100.0; - Real sigma_0 = 60.0; + Real C_T_ewp = interpolate_1d(wind_speed.dataPtr(), thrust_coeff.dataPtr(), z, wind_speed.size()); + + Real L_wake = std::pow(dx[0]*dx[1],0.5)/2.0; Real sigma_e = Vabs/(3.0*K_turb*L_wake)*(std::pow(2.0*K_turb*L_wake/Vabs + std::pow(sigma_0,2),3.0/2.0) - std::pow(sigma_0,3)); Real phi = std::atan2(v_vel(i,j,k),u_vel(i,j,k)); // Wind direction w.r.t the x-dreiction - Real fac = -std::pow(M_PI/8.0,0.5)*C_T_ewp*std::pow(R_ewp,2)*std::pow(Vabs,2)/(dx[0]*dx[1]*sigma_e)*std::exp(-0.5*std::pow((z - z_c_ewp)/sigma_e,2)); + Real fac = -std::pow(M_PI/8.0,0.5)*C_T_ewp*std::pow(rotor_rad,2)*std::pow(Vabs,2)/(dx[0]*dx[1]*sigma_e)*std::exp(-0.5*std::pow((z - hub_height)/sigma_e,2)); ewp_array(i,j,k,0) = fac*std::cos(phi)*Nturb_array(i,j,k); ewp_array(i,j,k,1) = fac*std::sin(phi)*Nturb_array(i,j,k); ewp_array(i,j,k,2) = 0.0; diff --git a/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp index 238c2de9f..1526a24f6 100644 --- a/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp +++ b/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp @@ -1,18 +1,19 @@ +#include +#include #include #include + using namespace amrex; -Real R = 30.0; -Real z_c = 100.0; -Real C_T = 0.8, C_TKE = 0.0; +Real C_TKE = 0.0; Real compute_A(Real z) { - Real d = std::min(std::fabs(z - z_c), R); - Real theta = std::acos(d/R); - Real A_s = R*R*theta - d*std::pow(R*R - d*d, 0.5); - Real A = M_PI*R*R/2.0 - A_s; + Real d = std::min(std::fabs(z - hub_height), rotor_rad); + Real theta = std::acos(d/rotor_rad); + Real A_s = rotor_rad*rotor_rad*theta - d*std::pow(rotor_rad*rotor_rad - d*d, 0.5); + Real A = M_PI*rotor_rad*rotor_rad/2.0 - A_s; return A; } @@ -23,7 +24,7 @@ Real compute_Aijk(Real z_k, Real z_kp1) Real A_k = compute_A(z_k); Real A_kp1 = compute_A(z_kp1); - Real check = (z_k - z_c)*(z_kp1 - z_c); + Real check = (z_k - hub_height)*(z_kp1 - hub_height); Real A_ijk; if(check > 0){ A_ijk = std::fabs(A_k -A_kp1); @@ -137,18 +138,20 @@ void fitch_source_terms_cellcentered (const Geometry& geom, // Compute Fitch source terms - Real Vabs = std::pow(u_vel(i,j,k)*u_vel(i,j,k) + - v_vel(i,j,k)*v_vel(i,j,k) + - w_vel(i,j,kk)*w_vel(i,j,kk), 0.5); + Real Vabs = std::pow(u_vel(i,j,k)*u_vel(i,j,k) + + v_vel(i,j,k)*v_vel(i,j,k) + + w_vel(i,j,kk)*w_vel(i,j,kk), 0.5); + + Real C_T = interpolate_1d(wind_speed.dataPtr(), thrust_coeff.dataPtr(), z, wind_speed.size()); - fitch_array(i,j,k,0) = Vabs; - fitch_array(i,j,k,1) = -0.5*Nturb_array(i,j,k)/(dx[0]*dx[1])*C_T*Vabs*Vabs*A_ijk/(z_kp1 - z_k); - fitch_array(i,j,k,2) = u_vel(i,j,k)/Vabs*fitch_array(i,j,k,1); - fitch_array(i,j,k,3) = v_vel(i,j,k)/Vabs*fitch_array(i,j,k,1); - fitch_array(i,j,k,4) = 0.5*Nturb_array(i,j,k)/(dx[0]*dx[1])*C_TKE*std::pow(Vabs,3)*A_ijk/(z_kp1 - z_k); + fitch_array(i,j,k,0) = Vabs; + fitch_array(i,j,k,1) = -0.5*Nturb_array(i,j,k)/(dx[0]*dx[1])*C_T*Vabs*Vabs*A_ijk/(z_kp1 - z_k); + fitch_array(i,j,k,2) = u_vel(i,j,k)/Vabs*fitch_array(i,j,k,1); + fitch_array(i,j,k,3) = v_vel(i,j,k)/Vabs*fitch_array(i,j,k,1); + fitch_array(i,j,k,4) = 0.5*Nturb_array(i,j,k)/(dx[0]*dx[1])*C_TKE*std::pow(Vabs,3)*A_ijk/(z_kp1 - z_k); //amrex::Gpu::Atomic::Add(sum_area, A_ijk); - }); + }); } //std::cout << "Checking sum here...." <<"\n"; //printf("%0.15g, %0.15g\n", *sum_area , M_PI*R*R); diff --git a/Source/WindFarmParametrization/WindFarm.H b/Source/WindFarmParametrization/WindFarm.H index 132e8706c..a0376ae52 100644 --- a/Source/WindFarmParametrization/WindFarm.H +++ b/Source/WindFarmParametrization/WindFarm.H @@ -5,6 +5,11 @@ #include #include +extern amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; +extern amrex::Vector wind_speed, thrust_coeff, power; + +void read_in_table(std::string windfarm_spec_table); + void advance_windfarm (int lev, const amrex::Geometry& geom, const amrex::Real& dt_advance, diff --git a/Source/WindFarmParametrization/WindFarm.cpp b/Source/WindFarmParametrization/WindFarm.cpp index 4162e7b73..3e54d02c3 100644 --- a/Source/WindFarmParametrization/WindFarm.cpp +++ b/Source/WindFarmParametrization/WindFarm.cpp @@ -3,6 +3,53 @@ #include using namespace amrex; +amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; +amrex::Vector wind_speed, thrust_coeff, power; + +void read_in_table(std::string windfarm_spec_table) +{ + + //The first line is the number of pairs entries for the power curve and thrust coefficient. + //The second line gives first the height in meters of the turbine hub, second, the diameter in + //meters of the rotor, third the standing thrust coefficient, and fourth the nominal power of + //the turbine in MW. + //The remaining lines contain the three values of: wind speed, thrust coefficient, and power production in kW. + + // Read turbine data from wind-turbine-1.tbl + std::ifstream file_turb_table(windfarm_spec_table); + if (!file_turb_table.is_open()) { + amrex::Error("Wind farm specifications table not found. Either the inputs is missing the " + "erf.windfarm_spec_table entry or the file specified in the entry - " + windfarm_spec_table + " is missing."); + } + else { + amrex::Print() << "Reading in wind farm specifications table: " << windfarm_spec_table << "\n"; + } + + int nlines; + file_turb_table >> nlines; + wind_speed.resize(nlines); + thrust_coeff.resize(nlines); + power.resize(nlines); + + Real rotor_dia; + file_turb_table >> hub_height >> rotor_dia >> thrust_coeff_standing >> nominal_power; + rotor_rad = rotor_dia*0.5; + if(rotor_rad > hub_height) { + amrex::Abort("The blade length is more than the hub height. Check the second line in wind-turbine-1.tbl. Aborting....."); + } + if(thrust_coeff_standing > 1.0) { + amrex::Abort("The standing thrust coefficient is greater than 1. Check the second line in wind-turbine-1.tbl. Aborting....."); + } + + for(int iline=0;iline> wind_speed[iline] >> thrust_coeff[iline] >> power[iline]; + if(thrust_coeff[iline] > 1.0) { + amrex::Abort("The thrust coefficient is greater than 1. Check wind-turbine-1.tbl. Aborting....."); + } + } + file_turb_table.close(); +} + void advance_windfarm (int lev, const Geometry& geom, diff --git a/Submodules/AMReX b/Submodules/AMReX index f141191fc..8255d7593 160000 --- a/Submodules/AMReX +++ b/Submodules/AMReX @@ -1 +1 @@ -Subproject commit f141191fc73cb22ee649f89569c1778fc910d1fb +Subproject commit 8255d75931e669766ac3748fd596daef76fb9d5e