diff --git a/examples/boussinesq/incompressible_eady.py b/examples/boussinesq/incompressible_eady.py new file mode 100644 index 000000000..4bf68f884 --- /dev/null +++ b/examples/boussinesq/incompressible_eady.py @@ -0,0 +1,233 @@ +""" +The Eady-Boussinesq problem in a vertical slice (one cell thick in the +y-direction) solved using the incompressible Boussinesq equations. The problem +is described in Yamazaki, Shipton, Cullen, Mitchell and Cotter, 2017: +``Vertical slice modelling of nonlinear Eady waves using a compatible finite +element method''. +""" + +from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter +from gusto import ( + Domain, IncompressibleEadyParameters, IncompressibleEadyEquations, + IncompressibleGeostrophicImbalance, SawyerEliassenU, Perturbation, + boussinesq_hydrostatic_balance, BoussinesqSolver, SemiImplicitQuasiNewton, + OutputParameters, IO, SSPRK3, DGUpwind, SUPGOptions, YComponent +) +from firedrake import ( + as_vector, SpatialCoordinate, solve, ds_t, ds_b, PeriodicRectangleMesh, + ExtrudedMesh, cos, sin, cosh, sinh, tanh, pi, Function, sqrt, TrialFunction, + TestFunction, inner, dx, div, FacetNormal, +) + +incompressible_eady_defaults = { + 'ncolumns': 300, + 'nlayers': 10, + 'dt': 300.0, + 'tmax': 11*24*60*60., + 'dumpfreq': 288, + 'dirname': 'incompressible_eady' +} + + +def incompressible_eady( + ncolumns=incompressible_eady_defaults['ncolumns'], + nlayers=incompressible_eady_defaults['nlayers'], + dt=incompressible_eady_defaults['dt'], + tmax=incompressible_eady_defaults['tmax'], + dumpfreq=incompressible_eady_defaults['dumpfreq'], + dirname=incompressible_eady_defaults['dirname'] +): + + # ------------------------------------------------------------------------ # + # Test case parameters + # ------------------------------------------------------------------------ # + + domain_width = 2000000. # Width of domain (m) + domain_thickness = 100000. # Thickness of domain (m) + domain_height = 10000. # Height of domain (m) + f = 1.e-04 # Coriolis parameter (1/s) + a = -4.5 # Amplitude of the perturbation (m/s) + Bu = 0.5 # Burger number + N = sqrt(2.5e-5) # Brunt-Vaisala frequency (1/s) + g = 10.0 # Acceleration due to gravity (m/s^2) + + # ------------------------------------------------------------------------ # + # Set up model objects + # ------------------------------------------------------------------------ # + + # Domain -- 2D periodic base mesh which is one cell thick + base_mesh = PeriodicRectangleMesh( + ncolumns, 1, domain_width, domain_thickness, quadrilateral=True) + mesh = ExtrudedMesh( + base_mesh, layers=nlayers, layer_height=domain_height/nlayers + ) + domain = Domain(mesh, dt, "RTCF", 1) + + # Equation + parameters = IncompressibleEadyParameters( + N=N, H=domain_height, L=0.5*domain_width, f=f, fourthorder=True, + deltax=domain_width/ncolumns, deltaz=domain_height/nlayers, Omega=0.5*f + ) + eqns = IncompressibleEadyEquations(domain, parameters) + + # I/O + output = OutputParameters(dirname=dirname, dumpfreq=dumpfreq) + diagnostic_fields = [ + YComponent('u'), IncompressibleGeostrophicImbalance(eqns), + SawyerEliassenU(eqns), Perturbation('p'), Perturbation('b') + ] + io = IO(domain, output, diagnostic_fields=diagnostic_fields) + + # Transport schemes and methods + b_opts = SUPGOptions() + transport_schemes = [ + SSPRK3(domain, "u"), + SSPRK3(domain, "b", options=b_opts) + ] + transport_methods = [ + DGUpwind(eqns, "u"), + DGUpwind(eqns, "b", ibp=b_opts.ibp) + ] + + # Linear solve + linear_solver = BoussinesqSolver(eqns) + + # Time stepper + stepper = SemiImplicitQuasiNewton( + eqns, io, transport_schemes, transport_methods, + linear_solver=linear_solver + ) + + # ------------------------------------------------------------------------ # + # Initial conditions + # ------------------------------------------------------------------------ # + + # Initial conditions + u0 = stepper.fields("u") + b0 = stepper.fields("b") + p0 = stepper.fields("p") + + # spaces + Vu = domain.spaces("HDiv") + Vb = domain.spaces("theta") + Vp = domain.spaces("DG") + + # parameters + x, _, z = SpatialCoordinate(mesh) + Nsq = parameters.Nsq + + # background buoyancy + bref = (z - 0.5*domain_height)*Nsq + b_b = Function(Vb).project(bref) + + # buoyancy perturbation + coth = lambda x: cosh(x)/sinh(x) + Z = lambda z: Bu*((z/domain_height)-0.5) + + n = Bu**(-1)*sqrt((Bu*0.5-tanh(Bu*0.5))*(coth(Bu*0.5)-Bu*0.5)) + + L = 0.5*domain_width + b_exp = a*sqrt(Nsq)*( + -(1.-Bu*0.5*coth(Bu*0.5))*sinh(Z(z))*cos(pi*(x-L)/L) + - n*Bu*cosh(Z(z))*sin(pi*(x-L)/L) + ) + b_pert = Function(Vb).interpolate(b_exp) + + # set total buoyancy + b0.project(b_b + b_pert) + + # calculate hydrostatic pressure + p_b = Function(Vp) + boussinesq_hydrostatic_balance(eqns, b_b, p_b) + boussinesq_hydrostatic_balance(eqns, b0, p0) + + # set x component of velocity + dbdy = parameters.dbdy + u = -dbdy/f*(z-0.5*domain_height) + + # set y component of velocity + v = Function(Vp).assign(0.) + + g = TrialFunction(Vu) + wg = TestFunction(Vu) + + n = FacetNormal(mesh) + + a = inner(wg, g)*dx + L = -div(wg)*p0*dx + inner(wg, n)*p0*(ds_t + ds_b) + pgrad = Function(Vu) + solve(a == L, pgrad) + + # get initial v + Vp = p0.function_space() + phi = TestFunction(Vp) + m = TrialFunction(Vp) + + a = f*phi*m*dx + L = phi*pgrad[0]*dx + solve(a == L, v) + + # set initial u + u_exp = as_vector([u, v, 0.]) + u0.project(u_exp) + + # set the background profiles + stepper.set_reference_profiles([('p', p_b), ('b', b_b)]) + + # ------------------------------------------------------------------------ # + # Run + # ------------------------------------------------------------------------ # + + stepper.run(t=0, tmax=tmax) + + +# ---------------------------------------------------------------------------- # +# MAIN +# ---------------------------------------------------------------------------- # + + +if __name__ == "__main__": + + parser = ArgumentParser( + description=__doc__, + formatter_class=ArgumentDefaultsHelpFormatter + ) + parser.add_argument( + '--ncolumns', + help="The number of columns in the vertical slice mesh.", + type=int, + default=incompressible_eady_defaults['ncolumns'] + ) + parser.add_argument( + '--nlayers', + help="The number of layers for the mesh.", + type=int, + default=incompressible_eady_defaults['nlayers'] + ) + parser.add_argument( + '--dt', + help="The time step in seconds.", + type=float, + default=incompressible_eady_defaults['dt'] + ) + parser.add_argument( + "--tmax", + help="The end time for the simulation in seconds.", + type=float, + default=incompressible_eady_defaults['tmax'] + ) + parser.add_argument( + '--dumpfreq', + help="The frequency at which to dump field output.", + type=int, + default=incompressible_eady_defaults['dumpfreq'] + ) + parser.add_argument( + '--dirname', + help="The name of the directory to write to.", + type=str, + default=incompressible_eady_defaults['dirname'] + ) + args, unknown = parser.parse_known_args() + + incompressible_eady(**vars(args)) diff --git a/examples/boussinesq/test_boussinesq_examples.py b/examples/boussinesq/test_boussinesq_examples.py index 53d3d479d..a4aac8d94 100644 --- a/examples/boussinesq/test_boussinesq_examples.py +++ b/examples/boussinesq/test_boussinesq_examples.py @@ -62,3 +62,16 @@ def test_skamarock_klemp_linear_bouss(): @pytest.mark.parallel(nprocs=2) def test_skamarock_klemp_linear_bouss_parallel(): test_skamarock_klemp_linear_bouss() + + +def test_incompressible_eady(): + from incompressible_eady import incompressible_eady + test_name = 'incompressible_eady' + incompressible_eady( + ncolumns=10, + nlayers=5, + dt=300.0, + tmax=1500.0, + dumpfreq=5, + dirname=make_dirname(test_name) + ) diff --git a/examples/compressible_euler/compressible_eady.py b/examples/compressible_euler/compressible_eady.py new file mode 100644 index 000000000..4b859a96e --- /dev/null +++ b/examples/compressible_euler/compressible_eady.py @@ -0,0 +1,243 @@ +""" +This solves the Eady problem in a vertical slice (one cell thick in the +y-direction) using the compressible Euler equations, following +Yamazaki and Cotter, 2025: +``A vertical slice frontogenesis test case for compressible nonhydrostatic +dynamical cores of atmospheric models''. +""" + +from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter +from gusto import ( + Domain, CompressibleEadyParameters, CompressibleEadyEquations, + Perturbation, compressible_hydrostatic_balance, + CompressibleSolver, SemiImplicitQuasiNewton, OutputParameters, IO, + SSPRK3, DGUpwind, SUPGOptions, YComponent, Exner +) +from gusto import thermodynamics as tde +from firedrake import ( + as_vector, SpatialCoordinate, solve, ds_b, ds_t, PeriodicRectangleMesh, + ExtrudedMesh, exp, cos, sin, cosh, sinh, tanh, pi, Function, sqrt, + TrialFunction, TestFunction, inner, dx, div, FacetNormal +) + +compressible_eady_defaults = { + 'ncolumns': 300, + 'nlayers': 10, + 'dt': 300.0, + 'tmax': 11*24*60*60., + 'dumpfreq': 288, + 'dirname': 'compressible_eady' +} + + +def compressible_eady( + ncolumns=compressible_eady_defaults['ncolumns'], + nlayers=compressible_eady_defaults['nlayers'], + dt=compressible_eady_defaults['dt'], + tmax=compressible_eady_defaults['tmax'], + dumpfreq=compressible_eady_defaults['dumpfreq'], + dirname=compressible_eady_defaults['dirname'] +): + + # ------------------------------------------------------------------------ # + # Test case parameters + # ------------------------------------------------------------------------ # + + domain_width = 2000000. # Width of domain (m) + domain_thickness = 100000. # Thickness of domain (m) + domain_height = 10000. # Height of domain (m) + f = 1.e-04 # Coriolis parameter (1/s) + a = -4.5 # Amplitude of the perturbation (m/s) + Bu = 0.5 # Burger number + N = sqrt(2.5e-5) # Brunt-Vaisala frequency (1/s) + g = 10.0 # Acceleration due to gravity (m/s^2) + Pi0 = 0.864 # Exner pressure at the surface + + # ------------------------------------------------------------------------ # + # Set up model objects + # ------------------------------------------------------------------------ # + + # Domain -- 2D periodic base mesh which is one cell thick + base_mesh = PeriodicRectangleMesh( + ncolumns, 1, domain_width, domain_thickness, quadrilateral=True) + mesh = ExtrudedMesh( + base_mesh, layers=nlayers, layer_height=domain_height/nlayers + ) + domain = Domain(mesh, dt, "RTCF", 1) + + # Equation + parameters = CompressibleEadyParameters( + N=N, H=domain_height, f=f, Pi0=Pi0, g=g, Omega=0.5*f + ) + eqns = CompressibleEadyEquations( + domain, parameters, u_transport_option='vector_advection_form' + ) + + # I/O + output = OutputParameters(dirname=dirname, dumpfreq=dumpfreq) + diagnostic_fields = [ + YComponent('u'), Exner(parameters), Perturbation('theta'), + Perturbation('Exner') + ] + io = IO(domain, output, diagnostic_fields=diagnostic_fields) + + # Transport schemes and methods + theta_opts = SUPGOptions() + transport_schemes = [ + SSPRK3(domain, "u"), + SSPRK3(domain, "rho"), + SSPRK3(domain, "theta", options=theta_opts) + ] + transport_methods = [ + DGUpwind(eqns, "u"), + DGUpwind(eqns, "rho"), + DGUpwind(eqns, "theta", ibp=theta_opts.ibp) + ] + + # Linear solver + linear_solver = CompressibleSolver(eqns) + + # Time stepper + stepper = SemiImplicitQuasiNewton( + eqns, io, transport_schemes, transport_methods, + linear_solver=linear_solver + ) + + # ------------------------------------------------------------------------ # + # Initial conditions + # ------------------------------------------------------------------------ # + + u0 = stepper.fields("u") + rho0 = stepper.fields("rho") + theta0 = stepper.fields("theta") + + # spaces + Vu = domain.spaces("HDiv") + Vt = domain.spaces("theta") + Vr = domain.spaces("DG") + + # first setup the background buoyancy profile + # z.grad(bref) = N**2 + x, _, z = SpatialCoordinate(mesh) + g = parameters.g + Nsq = parameters.Nsq + theta_surf = parameters.theta_surf + + # N^2 = (g/theta)dtheta/dz => dtheta/dz = theta N^2g => theta=theta_0exp(N^2gz) + theta_ref = theta_surf*exp(Nsq*(z - 0.5*domain_height)/g) + theta_b = Function(Vt).interpolate(theta_ref) + + # set theta_pert + coth = lambda x: cosh(x)/sinh(x) + Z = lambda z: Bu*((z/domain_height)-0.5) + + n = Bu**(-1)*sqrt((Bu*0.5-tanh(Bu*0.5))*(coth(Bu*0.5)-Bu*0.5)) + + L = 0.5*domain_width + theta_exp = a*theta_surf/g*sqrt(Nsq)*( + -(1.-Bu*0.5*coth(Bu*0.5))*sinh(Z(z))*cos(pi*(x-L)/L) + - n*Bu*cosh(Z(z))*sin(pi*(x-L)/L) + ) + theta_pert = Function(Vt).interpolate(theta_exp) + + # set theta0 + theta0.interpolate(theta_b + theta_pert) + + # calculate hydrostatic Pi + rho_b = Function(Vr) + compressible_hydrostatic_balance(eqns, theta_b, rho_b, solve_for_rho=True) + compressible_hydrostatic_balance(eqns, theta0, rho0, solve_for_rho=True) + + Pi = tde.exner_pressure(parameters, rho0, theta0) + + # set x component of velocity + cp = parameters.cp + dthetady = parameters.dthetady + u = cp*dthetady/f*(Pi-Pi0) + + # set y component of velocity by solving a problem + v = Function(Vr).assign(0.) + + # get Pi gradient + g = TrialFunction(Vu) + wg = TestFunction(Vu) + + n = FacetNormal(mesh) + + a = inner(wg, g)*dx + L = -div(wg)*Pi*dx + inner(wg, n)*Pi*(ds_t + ds_b) + pgrad = Function(Vu) + solve(a == L, pgrad) + + # get initial v + m = TrialFunction(Vr) + phi = TestFunction(Vr) + + a = phi*f*m*dx + L = phi*cp*theta0*pgrad[0]*dx + solve(a == L, v) + + # set initial u + u_exp = as_vector([u, v, 0.]) + u0.project(u_exp) + + # set the background profiles + stepper.set_reference_profiles([('rho', rho_b), ('theta', theta_b)]) + + # ------------------------------------------------------------------------ # + # Run + # ------------------------------------------------------------------------ # + + stepper.run(t=0, tmax=tmax) + + +# ---------------------------------------------------------------------------- # +# MAIN +# ---------------------------------------------------------------------------- # + + +if __name__ == "__main__": + + parser = ArgumentParser( + description=__doc__, + formatter_class=ArgumentDefaultsHelpFormatter + ) + parser.add_argument( + '--ncolumns', + help="The number of columns in the vertical slice mesh.", + type=int, + default=compressible_eady_defaults['ncolumns'] + ) + parser.add_argument( + '--nlayers', + help="The number of layers for the mesh.", + type=int, + default=compressible_eady_defaults['nlayers'] + ) + parser.add_argument( + '--dt', + help="The time step in seconds.", + type=float, + default=compressible_eady_defaults['dt'] + ) + parser.add_argument( + "--tmax", + help="The end time for the simulation in seconds.", + type=float, + default=compressible_eady_defaults['tmax'] + ) + parser.add_argument( + '--dumpfreq', + help="The frequency at which to dump field output.", + type=int, + default=compressible_eady_defaults['dumpfreq'] + ) + parser.add_argument( + '--dirname', + help="The name of the directory to write to.", + type=str, + default=compressible_eady_defaults['dirname'] + ) + args, unknown = parser.parse_known_args() + + compressible_eady(**vars(args)) diff --git a/examples/compressible_euler/test_compressible_euler_examples.py b/examples/compressible_euler/test_compressible_euler_examples.py index d74d5dade..c3a608275 100644 --- a/examples/compressible_euler/test_compressible_euler_examples.py +++ b/examples/compressible_euler/test_compressible_euler_examples.py @@ -121,3 +121,16 @@ def test_unsaturated_bubble(): @pytest.mark.parallel(nprocs=2) def test_unsaturated_bubble_parallel(): test_unsaturated_bubble() + + +def test_compressible_eady(): + from compressible_eady import compressible_eady + test_name = 'compressible_eady' + compressible_eady( + ncolumns=10, + nlayers=5, + dt=300.0, + tmax=1500.0, + dumpfreq=5, + dirname=make_dirname(test_name) + ) diff --git a/gusto/core/configuration.py b/gusto/core/configuration.py index 40852033a..a0290450e 100644 --- a/gusto/core/configuration.py +++ b/gusto/core/configuration.py @@ -7,8 +7,8 @@ __all__ = [ "Configuration", "IntegrateByParts", "TransportEquationType", "OutputParameters", - "BoussinesqParameters", "CompressibleParameters", - "ShallowWaterParameters", + "BoussinesqParameters", "CompressibleParameters", "ShallowWaterParameters", + "IncompressibleEadyParameters", "CompressibleEadyParameters", "EmbeddedDGOptions", "ConservativeEmbeddedDGOptions", "RecoveryOptions", "ConservativeRecoveryOptions", "SUPGOptions", "MixedFSOptions", "SpongeLayerParameters", "DiffusionParameters", "BoundaryLayerParameters", @@ -164,6 +164,38 @@ class ShallowWaterParameters(Configuration): q0 = None +class EadyParameters(Configuration): + """ + Base class of physical parameters for Eady problems + """ + dbdy = -1.0e-07 + H = None + L = None + f = None + deltax = None + deltaz = None + fourthorder = False + + +class IncompressibleEadyParameters(BoussinesqParameters, EadyParameters): + """ + Base class of physical parameters for incompressible Eady problems + """ + Nsq = BoussinesqParameters.N**2 + + +class CompressibleEadyParameters(CompressibleParameters, EadyParameters): + + """ + Physical parameters for Compressible Eady + """ + g = 10. + Nsq = CompressibleParameters.N**2 + theta_surf = 300. + dthetady = theta_surf/g*EadyParameters.dbdy + Pi0 = 0.0 + + class WrapperOptions(Configuration, metaclass=ABCMeta): """Base class for specifying options for a transport scheme.""" diff --git a/gusto/diagnostics/__init__.py b/gusto/diagnostics/__init__.py index 2d37128b9..337ceff39 100644 --- a/gusto/diagnostics/__init__.py +++ b/gusto/diagnostics/__init__.py @@ -1,3 +1,4 @@ from gusto.diagnostics.diagnostics import * # noqa +from gusto.diagnostics.boussinesq_diagnostics import * # noqa from gusto.diagnostics.shallow_water_diagnostics import * # noqa from gusto.diagnostics.compressible_euler_diagnostics import * # noqa \ No newline at end of file diff --git a/gusto/diagnostics/boussinesq_diagnostics.py b/gusto/diagnostics/boussinesq_diagnostics.py new file mode 100644 index 000000000..03e017e5c --- /dev/null +++ b/gusto/diagnostics/boussinesq_diagnostics.py @@ -0,0 +1,197 @@ +"""Some diagnostic fields for the Boussinesq equations.""" + +from firedrake import ( + as_vector, inner, dx, div, as_matrix, TrialFunction, TestFunction, + LinearVariationalProblem, LinearVariationalSolver, DirichletBC, Function, + grad, dot, FacetNormal, SpatialCoordinate, avg, jump, lhs, rhs, sqrt, + FunctionSpace, dS_v, dS_h +) +from gusto.diagnostics.diagnostics import DiagnosticField + +__all__ = [ + 'IncompressibleGeostrophicImbalance', 'SawyerEliassenU'] + + +class IncompressibleGeostrophicImbalance(DiagnosticField): + """Diagnostic for the amount of geostrophic imbalance.""" + name = "GeostrophicImbalance" + + def __init__(self, equations, space=None, method='interpolate'): + """ + Args: + equations (:class:`IncompressibleEadyEquations`): the equation set + being solved by the model. + space (:class:`FunctionSpace`, optional): the function space to + evaluate the diagnostic field in. Defaults to None, in which + case a default space will be chosen for this diagnostic. + method (str, optional): a string specifying the method of evaluation + for this diagnostic. Valid options are 'interpolate', 'project', + 'assign' and 'solve'. Defaults to 'interpolate'. + """ + self.equations = equations + self.parameters = equations.parameters + super().__init__(space=space, method=method, required_fields=("u", "b", "p")) + + def setup(self, domain, state_fields): + """ + Sets up the :class:`Function` for the diagnostic field + + Args: + domain (:class:`Domain`): the model's domain object. + state_fields (:class:`StateFields`): the model's field container. + """ + + u = state_fields("u") + b = state_fields("b") + p = state_fields("p") + f = self.parameters.f + Vu = u.function_space() + + v = TrialFunction(Vu) + w = TestFunction(Vu) + a = inner(w, v)*dx + L = (div(w)*p+inner(w, as_vector([f*u[1], 0.0, b])))*dx + + bcs = self.equations.bcs['u'] + + imbalance = Function(Vu) + self.expr = imbalance[0]/f + imbalanceproblem = LinearVariationalProblem(a, L, imbalance, bcs=bcs) + self.imbalance_solver = LinearVariationalSolver( + imbalanceproblem, solver_parameters={'ksp_type': 'cg'}) + super().setup(domain, state_fields) + + def compute(self): + """Compute the diagnostic field from the current state.""" + self.imbalance_solver.solve() + super().compute() + + +class SawyerEliassenU(DiagnosticField): + """ + Velocity associated with the Sawyer-Eliassen balance equation: the + secondary circulation associated with a stream function that ensures thermal + wind balance. + """ + name = "SawyerEliassenU" + + def __init__(self, equations): + """ + Args: + equations (:class:`IncompressibleEadyEquations`): the equation set + being solved by the model. + """ + space = equations.domain.spaces('HDiv') + self.parameters = equations.parameters + self.solve_implemented = True + super().__init__(space=space, method='solve', required_fields=("u", "b", "p")) + + def setup(self, domain, state_fields): + """ + Sets up the :class:`Function` for the diagnostic field + Args: + domain (:class:`Domain`): the model's domain object. + state_fields (:class:`StateFields`): the model's field container. + """ + + super().setup(domain, state_fields) + + u = state_fields("u") + b = state_fields("b") + v = inner(u, as_vector([0., 1., 0.])) + + # spaces + V0 = domain.spaces('H1') + Vu = domain.spaces('HDiv') + + # project b to V0 + b_v0 = Function(V0) + btri = TrialFunction(V0) + btes = TestFunction(V0) + a = inner(btes, btri) * dx + L = inner(btes, b) * dx + projectbproblem = LinearVariationalProblem(a, L, b_v0) + self.project_b_solver = LinearVariationalSolver( + projectbproblem, solver_parameters={'ksp_type': 'cg'}) + + # project v to V0 + v_v0 = Function(V0) + vtri = TrialFunction(V0) + vtes = TestFunction(V0) + a = inner(vtes, vtri) * dx + L = inner(vtes, v) * dx + projectvproblem = LinearVariationalProblem(a, L, v_v0) + self.project_v_solver = LinearVariationalSolver( + projectvproblem, solver_parameters={'ksp_type': 'cg'}) + + # stm/psi is a stream function + stm = Function(V0) + psi = TrialFunction(V0) + xsi = TestFunction(V0) + + f = self.parameters.f + H = self.parameters.H + L = self.parameters.L + dbdy = self.parameters.dbdy + _, _, z = SpatialCoordinate(domain.mesh) + + bcs = [DirichletBC(V0, 0., "bottom"), + DirichletBC(V0, 0., "top")] + + Mat = as_matrix([[b.dx(2), 0., -f*v_v0.dx(2)], + [0., 0., 0.], + [-b_v0.dx(0), 0., f**2+f*v_v0.dx(0)]]) + + Equ = ( + inner(grad(xsi), dot(Mat, grad(psi))) + - dbdy*inner(grad(xsi), as_vector([-v, 0., f*(z-H/2)])) + )*dx + + # fourth-order terms + if self.parameters.fourthorder: + R = FunctionSpace(domain.mesh, "R", 0) + eps = Function(R, val=0.0001) + brennersigma = Function(R, val=10.0) + n = FacetNormal(domain.mesh) + deltax = self.parameters.deltax + deltaz = self.parameters.deltaz + + nn = as_matrix([[sqrt(brennersigma/deltax), 0., 0.], + [0., 0., 0.], + [0., 0., sqrt(brennersigma/deltaz)]]) + + mu = as_matrix([[1., 0., 0.], + [0., 0., 0.], + [0., 0., H/L]]) + + # anisotropic form + Equ += eps*( + div(dot(mu, grad(psi)))*div(dot(mu, grad(xsi)))*dx + - ( + avg(dot(dot(grad(grad(psi)), n), n))*jump(grad(xsi), n=n) + + avg(dot(dot(grad(grad(xsi)), n), n))*jump(grad(psi), n=n) + - jump(nn*grad(psi), n=n)*jump(nn*grad(xsi), n=n) + )*(dS_h + dS_v) + ) + + Au = lhs(Equ) + Lu = rhs(Equ) + stmproblem = LinearVariationalProblem(Au, Lu, stm, bcs=bcs) + self.stream_function_solver = LinearVariationalSolver( + stmproblem, solver_parameters={'ksp_type': 'cg'}) + + # solve for sawyer_eliassen u + utrial = TrialFunction(Vu) + w = TestFunction(Vu) + a = inner(w, utrial)*dx + L = (w[0]*(-stm.dx(2))+w[2]*(stm.dx(0)))*dx + ugproblem = LinearVariationalProblem(a, L, self.field) + self.sawyer_eliassen_u_solver = LinearVariationalSolver( + ugproblem, solver_parameters={'ksp_type': 'cg'}) + + def compute(self): + """Compute the diagnostic field from the current state.""" + self.project_b_solver.solve() + self.project_v_solver.solve() + self.stream_function_solver.solve() + self.sawyer_eliassen_u_solver.solve() diff --git a/gusto/equations/boussinesq_equations.py b/gusto/equations/boussinesq_equations.py index b9823f110..3f3cc3ee9 100644 --- a/gusto/equations/boussinesq_equations.py +++ b/gusto/equations/boussinesq_equations.py @@ -1,6 +1,9 @@ """Defines the Boussinesq equations.""" -from firedrake import inner, dx, div, cross, split, as_vector +from firedrake import ( + inner, dx, div, cross, split, as_vector, Function, TestFunctions, + SpatialCoordinate +) from firedrake.fml import subject from gusto.core.labels import ( time_derivative, transport, prognostic, linearisation, @@ -13,7 +16,10 @@ ) from gusto.equations.prognostic_equations import PrognosticEquationSet -__all__ = ["BoussinesqEquations", "LinearBoussinesqEquations"] +__all__ = [ + "BoussinesqEquations", "LinearBoussinesqEquations", + "IncompressibleEadyEquations" +] class BoussinesqEquations(PrognosticEquationSet): @@ -284,3 +290,72 @@ def __init__(self, domain, parameters, # Use the underlying routine to do a first linearisation of # the equations self.linearise_equation_set() + + +class IncompressibleEadyEquations(BoussinesqEquations): + """ + The incompressible Boussinesq equations in a vertical slice, augmented with + a third velocity component normal to the slice, as used in the Eady problem. + """ + + def __init__(self, domain, parameters, + space_names=None, linearisation_map='default', + u_transport_option="vector_invariant_form", + no_normal_flow_bc_ids=None): + """ + Args: + domain (:class:`Domain`): the model's domain object, containing the + mesh and the compatible function spaces. + parameters (:class:`Configuration`, optional): an object containing + the model's physical parameters. + space_names (dict, optional): a dictionary of strings for names of + the function spaces to use for the spatial discretisation. The + keys are the names of the prognostic variables. Defaults to None + in which case the spaces are taken from the de Rham complex. + linearisation_map (func, optional): a function specifying which + terms in the equation set to linearise. If None is specified + then no terms are linearised. Defaults to the string 'default', + in which case the linearisation includes time derivatives and + scalar transport terms. + u_transport_option (str, optional): specifies the transport term + used for the velocity equation. Supported options are: + 'vector_invariant_form', 'vector_advection_form' and + 'circulation_form'. + Defaults to 'vector_invariant_form'. + no_normal_flow_bc_ids (list, optional): a list of IDs of domain + boundaries at which no normal flow will be enforced. Defaults to + None. + """ + + if linearisation_map == 'default': + # Default linearisation is time derivatives, pressure gradient, + # Coriolis and transport term from depth equation + linearisation_map = lambda t: \ + (any(t.has_label(time_derivative, pressure_gradient, coriolis, + gravity, divergence, incompressible)) + or (t.get(prognostic) in ['p', 'b'] and t.has_label(transport))) + super().__init__(domain=domain, + parameters=parameters, + compressible=False, + space_names=space_names, + linearisation_map=linearisation_map, + u_transport_option=u_transport_option, + no_normal_flow_bc_ids=no_normal_flow_bc_ids, + active_tracers=None) + + dbdy = parameters.dbdy + H = parameters.H + _, _, z = SpatialCoordinate(domain.mesh) + eady_exp = Function(domain.spaces("DG")).interpolate(z-H/2.) + y_vec = as_vector([0., 1., 0.]) + + W = self.function_space + w, _, gamma = TestFunctions(W) + X = self.X + u, _, _ = split(X) + + self.residual += subject(prognostic( + dbdy*eady_exp*inner(w, y_vec)*dx, "u"), X) + + self.residual += subject(prognostic( + gamma*dbdy*inner(u, y_vec)*dx, "b"), X) diff --git a/gusto/equations/compressible_euler_equations.py b/gusto/equations/compressible_euler_equations.py index 73bc76d45..73043d2ca 100644 --- a/gusto/equations/compressible_euler_equations.py +++ b/gusto/equations/compressible_euler_equations.py @@ -2,7 +2,7 @@ from firedrake import ( sin, pi, inner, dx, div, cross, FunctionSpace, FacetNormal, jump, avg, dS_v, - conditional, SpatialCoordinate, split, Constant, as_vector + conditional, SpatialCoordinate, split, Constant, as_vector, TestFunctions ) from firedrake.fml import subject, replace_subject from gusto.core.labels import ( @@ -18,7 +18,10 @@ from gusto.equations.active_tracers import Phases, TracerVariableType from gusto.equations.prognostic_equations import PrognosticEquationSet -__all__ = ["CompressibleEulerEquations", "HydrostaticCompressibleEulerEquations"] +__all__ = [ + "CompressibleEulerEquations", "HydrostaticCompressibleEulerEquations", + "CompressibleEadyEquations" +] class CompressibleEulerEquations(PrognosticEquationSet): @@ -382,3 +385,72 @@ def hydrostatic_projection(self, term, field_name): # to replace the subject with the new hydrostatic subject # - add the hydrostatic label return replace_subject(new_subj, old_idx=f_idx)(term) + + +class CompressibleEadyEquations(CompressibleEulerEquations): + """ + The compressible equations in a vertical slice, augmented with a third + velocity component normal to the slice, as used in the Eady problem. + """ + + def __init__(self, domain, parameters, sponge_options=None, + extra_terms=None, space_names=None, + linearisation_map='default', + u_transport_option="vector_invariant_form", + no_normal_flow_bc_ids=None): + """ + Args: + domain (:class:`Domain`): the model's domain object, containing the + mesh and the compatible function spaces. + parameters (:class:`Configuration`, optional): an object containing + the model's physical parameters. + sponge_options (:class:`SpongeLayerParameters`, optional): any + parameters for applying a sponge layer to the upper boundary. + Defaults to None. + extra_terms (:class:`ufl.Expr`, optional): any extra terms to be + included in the equation set. Defaults to None. + space_names (dict, optional): a dictionary of strings for names of + the function spaces to use for the spatial discretisation. The + keys are the names of the prognostic variables. Defaults to None + in which case the spaces are taken from the de Rham complex. + linearisation_map (func, optional): a function specifying which + terms in the equation set to linearise. If None is specified + then no terms are linearised. Defaults to the string 'default', + in which case the linearisation includes time derivatives and + scalar transport terms. + u_transport_option (str, optional): specifies the transport term + used for the velocity equation. Supported options are: + 'vector_invariant_form', 'vector_advection_form' and + 'circulation_form'. + Defaults to 'vector_invariant_form'. + no_normal_flow_bc_ids (list, optional): a list of IDs of domain + boundaries at which no normal flow will be enforced. Defaults to + None. + """ + + super().__init__(domain, parameters, + sponge_options=sponge_options, + extra_terms=extra_terms, space_names=space_names, + linearisation_map=linearisation_map, + u_transport_option=u_transport_option, + diffusion_options=None, + no_normal_flow_bc_ids=no_normal_flow_bc_ids, + active_tracers=None) + + dthetady = parameters.dthetady + cp = parameters.cp + y_vec = as_vector([0., 1., 0.]) + + W = self.function_space + w, _, gamma = TestFunctions(W) + X = self.X + u, rho, theta = split(X) + + exner_pi = exner_pressure(parameters, rho, theta) + Pi0 = self.parameters.Pi0 + + self.residual -= subject(prognostic( + cp*dthetady*(exner_pi - Pi0)*inner(w, y_vec)*dx, "u"), X) + + self.residual += subject(prognostic( + gamma*(dthetady*inner(u, y_vec))*dx, "theta"), X) diff --git a/gusto/solvers/linear_solvers.py b/gusto/solvers/linear_solvers.py index 213d09df2..c05bc3c99 100644 --- a/gusto/solvers/linear_solvers.py +++ b/gusto/solvers/linear_solvers.py @@ -515,7 +515,7 @@ def V(u): eqn += dt*self.equations.mu*inner(w, k)*inner(u, k)*dx if equation.parameters.Omega is not None: - Omega = as_vector((0, 0, equation.parameter.Omega)) + Omega = as_vector((0, 0, equation.parameters.Omega)) eqn += inner(w, cross(2*Omega, u))*dx aeqn = lhs(eqn)