Integration
PinCFlow.Integration
— ModuleIntegration
Module for integration of the full system.
Provides helper functions for computing the time step, managing time levels and synchronizing fields, as well as the main function for running PinCFlow.jl.
See also
PinCFlow.Integration.backup_predictands
— Functionbackup_predictands(state::State)::Tuple{<:Predictands, <:TracerPredictands}
Return a tuple with a copy of the predictands and tracer predictands.
Arguments
state
: Model state
PinCFlow.Integration.compute_time_step
— Functioncompute_time_step(state::State)::AbstractFloat
Compute and return an adaptive time step based on several stability criteria.
If state.namelists.discretization.adaptive_time_step
is set to true
, the returned time step is given by
\[\Delta t = \min \left(\Delta t_\mathrm{CFL}, \Delta t_\mathrm{WKB}, \Delta t_\mathrm{viscous}, \Delta t_{\max}\right),\]
where $\Delta t_\mathrm{CFL}$ and $\Delta t_\mathrm{WKB}$ are computed from CFL conditions with respect to the resolved flow and unresolved gravity waves, respectively, $\Delta t_\mathrm{viscous}$ is determined from a von Neumann condition that takes the viscosity into account and $\Delta t_{\max}$ is an upper limit specified in state.namelists.discretization
. Otherwise, the returned time step is equal to $\Delta t_{\max}$. If $\Delta t$ is smaller than $\Delta t_{\min}$ (also specified in state.namelists.discretization
), an error is thrown.
The individual stability criteria are as follows.
CFL condition with respect to the resolved flow (where $w$ is computed with
compute_vertical_wind
):\[\Delta t_\mathrm{CFL} = \mu_\mathrm{CFL} \min\limits_\mathrm{global} \left[\frac{\Delta \widehat{x}}{u_{\max}}, \frac{\Delta \widehat{y}}{v_{\max}}, \min \left(\frac{J \Delta \widehat{z}}{w}\right)\right]\]
CFL condition with respect to the group velocities of unresolved gravity waves (where $J_{\min}$ is the minimum Jacobian in a one-grid-cell radius and $c_{\mathrm{g} z}$ is the maximum vertical group velocity within a grid cell):
\[\Delta t_\mathrm{WKB} = \mu_\mathrm{WKB} \min\limits_\mathrm{global} \left[\frac{\Delta \widehat{x}}{c_{\mathrm{g} x, \max}}, \frac{\Delta \widehat{y}}{c_{\mathrm{g} y, \max}}, \min \left(\frac{J_{\min} \Delta \widehat{z}}{c_{\mathrm{g} z}}\right)\right]\]
Von Neumann condition (with $\mathrm{Re}$ being the Reynolds number):
\[\Delta t_\mathrm{viscous} = \frac{\mathrm{Re}}{2} \min\limits_\mathrm{global} \left[\left(\Delta \widehat{x}\right)^2, \left(\Delta \widehat{y}\right)^2, \left(J \Delta \widehat{z}\right)^2\right]\]
Arguments
state
: Model state.
See also
PinCFlow.Integration.explicit_integration!
— Functionexplicit_integration!(
state::State,
p0::Predictands,
dtstage::AbstractFloat,
time::AbstractFloat,
side::LHS,
)
Integrate the left-hand sides of the prognostic equations with a Runge-Kutta time step.
At each Runge-Kutta stage, the prognostic variables are first reconstructed and their advective and diffusive fluxes are calculated. Subsequently, each variable is updated with its integrated left-hand side, followed immediately by an implicit Euler step (the size of which is the fractional time step at the current Runge-Kutta stage) that accounts for the Rayleigh-damping imposed by the LHS sponge. After the Runge-Kutta loop, if the atmosphere is compressible, the Exner-pressure is updated with a full implicit Euler step to account for it being impacted by the Rayleigh damping of the mass-weighted potential temperature.
explicit_integration!(
state::State,
p0::Predictands,
dtstage::AbstractFloat,
time::AbstractFloat,
side::RHS,
)
Perform an explicit Euler step on the right-hand sides of the prognostic equations and, if the atmosphere is compressible, update the Exner-pressure such that it is synchronized with the mass-weighted potential temperature.
Arguments
state
: Model state.p0
: The predictands that are used to compute the transporting velocities in the computation of the fluxes.dtstage
: Fractional time step.time
: Simulation time.side
: Side of the equations.
PinCFlow.Integration.implicit_integration!
— Functionimplicit_integration!(
state::State,
dtstage::AbstractFloat,
time::AbstractFloat,
ntotalbicg::Integer,
side::RHS,
rayleigh_factor::AbstractFloat,
iout::Integer,
machine_start_time::DateTime,
)
Perform an implicit Euler step on the right-hand sides of the prognostic equations, solve the Poisson equation and correct the Exner-pressure, momentum and density fluctuations accordingly.
Arguments
state
: Model state.dtstage
: Fractional time step.time
: Simulation time.ntotalbicg
: BicGStab-iterations counter.side
: Side of the equations.rayleigh_factor
: Factor by which the Rayleigh-damping coefficient is multiplied.iout
: Output counter.machine_start_time
: Wall-clock start time.
PinCFlow.Integration.integrate
— Functionintegrate(namelists::Namelists)
Initialize the model state and integrate it in time.
This method performs the complete time integration of the governing equations, using a semi-implicit time stepping scheme. It handles initialization, time stepping and output of the simulation data.
The initialization begins with the construction of the model state (an instance of the composite type State
), which involves the setup of the MPI parallelization and the definition of all arrays that are needed repeatedly during the simulation. This is followed by an (optional) initial cleaning, in which the Poisson solver is called to ensure that the initial dynamic fields satisfy the divergence constraint imposed by the thermodynamic energy equation. Afterwards, the initialization of MSGWaM is completed by adding ray volumes to the previously defined arrays. If the simulation is supposed to start from a previous model state, the fields are then overwritten with the data in the corresponding input file. Finally, the output file is created and the initial state is written into it.
At the beginning of each time-loop iteration, the time step is determined from several stability criteria, using compute_time_step
. In case the updated simulation time is later than the next output time, the time step is corrected accordingly. Subsequently, the damping coefficients of the sponges (which may depend on the time step) are calculated. Following this, MSGWaM updates the unresolved gravity-wave field and computes the corresponding mean-flow impact. Afterwards, the resolved flow is updated in a semi-implicit time step, comprised of the following stages.
Explicit RK3 integration of LHS over $\Delta t / 2$.
Implicit Euler integration of RHS over $\Delta t / 2$.
Explicit Euler integration of RHS over $\Delta t / 2$.
Explicit RK3 integration of LHS over $\Delta t$.
Implicit Euler integration of RHS over $\Delta t / 2$.
Therein, the left-hand sides of the equations include advective fluxes, diffusion terms, rotation and heating, whereas the pressure gradient, buoyancy term and momentum-flux divergence due to unresolved gravity waves are on the right-hand sides. Boundary conditions are enforced continuously. At the end of the time step, the updated fields are written into the output file if the next output time has been reached.
Arguments
namelists
: Namelists with all model parameters.
See also
PinCFlow.Integration.modify_compressible_wind!
— Functionmodify_compressible_wind!(state::State, operation::Function)
Modify the wind with $J P$ if the atmosphere is compressible by dispatching to the appropriate method.
modify_compressible_wind!(
state::State,
operation::Function,
model::Union{Boussinesq, PseudoIncompressible},
)
Return in non-compressible modes.
modify_compressible_wind!(
state::State,
operation::Function,
model::Compressible,
)
Interpolate $J P$ to the wind grids and replace the wind components with the result of applying operation
to them and the interpolations.
Arguments
state
: Model state.operation
: Binary operation used for modification.model
: Dynamic equations.
PinCFlow.Integration.reset_predictands!
— Functionreset_predictands!(
state::State,
predictands::Predictands,
tracerpredictands::TracerPredictands,
)
Reset fields in state
to those in predictands
and tracerpredictands
by dispatching to specific methods.
reset_predictands!(state::State, tracerpredictands::TracerPredictands)
Reset fields in state.tracer.tracerpredictands
to those in tracerpredictands
.
reset_predictands!(
state::State,
predictands::Predictands,
model::Union{Boussinesq, PseudoIncompressible},
)
Reset the density, density fluctuations and wind components in state.variables.predictands
to those in predictands
.
reset_predictands!(state::State, predictands::Predictands, model::Compressible)
Reset the density, density fluctuations, wind components, Exner pressure and mass-weighted potential temperature (i.e. all fields) in state.variables.predictands
to those in predictands
.
Arguments
state
: Model state.predictands
: Fields to reset to.tracerpredictands
: Tracer fields to reset to.model
: Dynamic equations.
PinCFlow.Integration.save_backups!
— Functionsave_backups!(state::State, variables::Vararg{Symbol})
Copy the specified fields in state.variables.predictands
to their counterparts in state.variables.backups
.
Arguments
state
: Model state.variables
: Names of the fields to create backups of.
PinCFlow.Integration.synchronize_compressible_atmosphere!
— Functionsynchronize_compressible_atmosphere!(state::State, predictands::Predictands)
Synchronize state.atmosphere.pbar
with predictands.p
if the atmosphere is compressible by dispatching to the appropriate method.
synchronize_compressible_atmosphere!(
state::State,
predictands::Predictands,
model::Union{Boussinesq, PseudoIncompressible},
)
Return in non-compressible modes.
synchronize_compressible_atmosphere!(
state::State,
predictands::Predictands,
model::Compressible,
)
Synchronize state.atmosphere.pbar
with predictands.p
.
In compressible mode, $P$ is time-dependent. In the update of $P$, only state.variables.predictands.p
is changed, so that the old values of $P$ are retained in the respective field of state.atmosphere
. When these are no longer needed, this method is used to update the field accordingly.
Arguments
state
: Model state.predictands
: Predictands to use for the synchronization of the mass-weighted potential temperature.model
: Dynamic equations.
PinCFlow.Integration.synchronize_density_fluctuations!
— Functionsynchronize_density_fluctuations!(state::State)
Synchronize the density fluctuations in state.variables.predictands.rhop
with the density in state.variables.predictands.rho
by dispatching to a model-specific method.
synchronize_density_fluctuations!(state::State, model::Boussinesq)
Return in Boussinesq mode.
In Boussinesq mode, density fluctuations don't require synchronization, since the density is assumed constant except in the buoyancy equation.
synchronize_density_fluctuations!(state::State, model::PseudoIncompressible)
Synchronize the density fluctuations in state.variables.predictands.rhop
with the density in state.variables.predictands.rho
.
The density fluctuations are defined as the product of the mass-weighted potential temperature and the fluctuations of the inverse potential temperature. In pseudo-incompressible mode, $P$ is constant, so that this is reduced to the difference between $\rho$ and $\overline{\rho}$, i.e.
\[\rho' = \frac{P}{\theta} - \frac{P}{\overline{\theta}} = \rho - \overline{\rho}.\]
synchronize_density_fluctuations!(state::State, model::Compressible)
Synchronize the density fluctuations in state.variables.predictands.rhop
with the density in state.variables.predictands.rho
.
In compressible mode, $P$ is time-dependent, so that the density fluctuations are not reduced to the difference between $\rho$ and $\overline{\rho}$, i.e.
\[\rho' = \frac{P}{\theta} - \frac{P}{\overline{\theta}} = \rho - \frac{P}{\overline{\theta}}.\]
Arguments
state
: Model state.model
: Dynamic equations.
PinCFlow.Integration.wkb_integration!
— Functionwkb_integration!(state::State, dtstage::AbstractFloat)
Use MSGWaM to update the properties of the unresolved gravity-wave field and compute its impact on the resolved flow.
In the first step, MSGWaM's saturation scheme is applied to account for the impact of wave breaking on the wave-action field. Next, the ray volumes are propagated via integration of the ray equations with a Runge-Kutta time step. Ray volumes that have grown larger than the cells of the model grid are then split before their array positions are shifted such that they are attributed to the correct grid cells. Afterwards, ray volumes are merged in cells where their count exceeds a specified threshold. Finally, the ray volumes in boundary and halo cells are updated and the mean-flow impact is calculated.
Arguments
state
: Model state.dtstage
: Time step.