diff --git a/examples/compressible_euler/mountain_hydrostatic.py b/examples/compressible_euler/schaer_mountain.py similarity index 62% rename from examples/compressible_euler/mountain_hydrostatic.py rename to examples/compressible_euler/schaer_mountain.py index 1fc113a71..d60679d1d 100644 --- a/examples/compressible_euler/mountain_hydrostatic.py +++ b/examples/compressible_euler/schaer_mountain.py @@ -1,68 +1,71 @@ """ -The hydrostatic 1 metre high mountain test case from Melvin et al, 2010: -``An inherently mass-conserving iterative semi-implicit semi-Lagrangian -discretization of the non-hydrostatic vertical-slice equations.'', QJRMS. +The Schär mountain test case of Schär et al, 2002: +``A new terrain-following vertical coordinate formulation for atmospheric +prediction models.'', MWR. -This test describes a wave over a mountain in a hydrostatic atmosphere. +This test describes a wave over a set of idealised mountains, testing how the +discretisation handles orography. The setup used here uses the order 1 finite elements. """ - from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter + from firedrake import ( as_vector, VectorFunctionSpace, PeriodicIntervalMesh, ExtrudedMesh, - SpatialCoordinate, exp, pi, cos, Function, conditional, Mesh, Constant + SpatialCoordinate, exp, pi, cos, Function, Mesh, Constant ) from gusto import ( - Domain, IO, OutputParameters, SemiImplicitQuasiNewton, SSPRK3, DGUpwind, - TrapeziumRule, SUPGOptions, ZComponent, Perturbation, - CompressibleParameters, HydrostaticCompressibleEulerEquations, - CompressibleSolver, compressible_hydrostatic_balance, HydrostaticImbalance, - SpongeLayerParameters, MinKernel, MaxKernel, logger + Domain, CompressibleParameters, CompressibleSolver, logger, + OutputParameters, IO, SSPRK3, DGUpwind, SemiImplicitQuasiNewton, + compressible_hydrostatic_balance, SpongeLayerParameters, Exner, ZComponent, + Perturbation, SUPGOptions, TrapeziumRule, MaxKernel, MinKernel, + CompressibleEulerEquations, SubcyclingOptions, RungeKuttaFormulation ) -mountain_hydrostatic_defaults = { - 'ncolumns': 200, - 'nlayers': 120, - 'dt': 5.0, - 'tmax': 15000., - 'dumpfreq': 1500, - 'dirname': 'mountain_hydrostatic' +schaer_mountain_defaults = { + 'ncolumns': 100, + 'nlayers': 50, + 'dt': 8.0, + 'tmax': 5*60*60., # 5 hours + 'dumpfreq': 2250, # dump at end with default settings + 'dirname': 'schaer_mountain' } -def mountain_hydrostatic( - ncolumns=mountain_hydrostatic_defaults['ncolumns'], - nlayers=mountain_hydrostatic_defaults['nlayers'], - dt=mountain_hydrostatic_defaults['dt'], - tmax=mountain_hydrostatic_defaults['tmax'], - dumpfreq=mountain_hydrostatic_defaults['dumpfreq'], - dirname=mountain_hydrostatic_defaults['dirname'] +def schaer_mountain( + ncolumns=schaer_mountain_defaults['ncolumns'], + nlayers=schaer_mountain_defaults['nlayers'], + dt=schaer_mountain_defaults['dt'], + tmax=schaer_mountain_defaults['tmax'], + dumpfreq=schaer_mountain_defaults['dumpfreq'], + dirname=schaer_mountain_defaults['dirname'] ): # ------------------------------------------------------------------------ # # Parameters for test case # ------------------------------------------------------------------------ # - domain_width = 240000. # width of domain in x direction, in m - domain_height = 50000. # height of model top, in m - a = 10000. # scale width of mountain, in m - hm = 1. # height of mountain, in m - zh = 5000. # height at which mesh is no longer distorted, in m - Tsurf = 250. # temperature of surface, in K - initial_wind = 20.0 # initial horizontal wind, in m/s - sponge_depth = 20000.0 # depth of sponge layer, in m - g = 9.80665 # acceleration due to gravity, in m/s^2 - cp = 1004. # specific heat capacity at constant pressure - sponge_mu = 0.15 # parameter for strength of sponge layer, in J/kg/K + domain_width = 100000. # width of domain in x direction, in m + domain_height = 30000. # height of model top, in m + a = 5000. # scale width of mountain profile, in m + lamda = 4000. # scale width of individual mountains, in m + hm = 250. # height of mountain, in m + Tsurf = 288. # temperature of surface, in K + initial_wind = 10.0 # initial horizontal wind, in m/s + sponge_depth = 10000.0 # depth of sponge layer, in m + g = 9.810616 # acceleration due to gravity, in m/s^2 + cp = 1004.5 # specific heat capacity at constant pressure + mu_dt = 1.2 # strength of sponge layer, no units exner_surf = 1.0 # maximum value of Exner pressure at surface - max_iterations = 10 # maximum number of hydrostatic balance iterations - tolerance = 1e-7 # tolerance for hydrostatic balance iteration + max_iterations = 20 # maximum number of hydrostatic balance iterations + tolerance = 1e-8 # tolerance for hydrostatic balance iteration # ------------------------------------------------------------------------ # # Our settings for this set up # ------------------------------------------------------------------------ # + spinup_steps = 5 # Not necessary but helps balance initial conditions + alpha = 0.51 # Necessary to absorb grid scale waves element_order = 1 u_eqn_type = 'vector_invariant_form' @@ -81,9 +84,9 @@ def mountain_hydrostatic( # Describe the mountain xc = domain_width/2. x, z = SpatialCoordinate(ext_mesh) - zs = hm * a**2 / ((x - xc)**2 + a**2) + zs = hm * exp(-((x - xc)/a)**2) * (cos(pi*(x - xc)/lamda))**2 xexpr = as_vector( - [x, conditional(z < zh, z + cos(0.5 * pi * z / zh)**6 * zs, z)] + [x, z + ((domain_height - z) / domain_height) * zs] ) # Make new mesh @@ -95,64 +98,50 @@ def mountain_hydrostatic( # Equation parameters = CompressibleParameters(g=g, cp=cp) sponge = SpongeLayerParameters( - H=domain_height, z_level=domain_height-sponge_depth, mubar=sponge_mu/dt + H=domain_height, z_level=domain_height-sponge_depth, mubar=mu_dt/dt ) - eqns = HydrostaticCompressibleEulerEquations( + eqns = CompressibleEulerEquations( domain, parameters, sponge_options=sponge, u_transport_option=u_eqn_type ) # I/O output = OutputParameters( - dirname=dirname, dumpfreq=dumpfreq, dump_vtus=True, dump_nc=False + dirname=dirname, dumpfreq=dumpfreq, dump_vtus=False, dump_nc=True ) diagnostic_fields = [ - ZComponent('u'), HydrostaticImbalance(eqns), - Perturbation('theta'), Perturbation('rho') + Exner(parameters), ZComponent('u'), Perturbation('theta'), + Perturbation('rho') ] io = IO(domain, output, diagnostic_fields=diagnostic_fields) # Transport schemes + subcycling_opts = SubcyclingOptions(subcycle_by_courant=0.25) theta_opts = SUPGOptions() transported_fields = [ - TrapeziumRule(domain, "u"), - SSPRK3(domain, "rho"), - SSPRK3(domain, "theta", options=theta_opts) + TrapeziumRule(domain, "u", subcycling_options=subcycling_opts), + SSPRK3( + domain, "rho", rk_formulation=RungeKuttaFormulation.predictor, + subcycling_options=subcycling_opts + ), + SSPRK3( + domain, "theta", options=theta_opts, + subcycling_options=subcycling_opts + ) ] transport_methods = [ DGUpwind(eqns, "u"), - DGUpwind(eqns, "rho"), + DGUpwind(eqns, "rho", advective_then_flux=True), DGUpwind(eqns, "theta", ibp=theta_opts.ibp) ] # Linear solver - params = {'mat_type': 'matfree', - 'ksp_type': 'preonly', - 'pc_type': 'python', - 'pc_python_type': 'firedrake.SCPC', - # Velocity mass operator is singular in the hydrostatic case. - # So for reconstruction, we eliminate rho into u - 'pc_sc_eliminate_fields': '1, 0', - 'condensed_field': {'ksp_type': 'fgmres', - 'ksp_rtol': 1.0e-8, - 'ksp_atol': 1.0e-8, - 'ksp_max_it': 100, - 'pc_type': 'gamg', - 'pc_gamg_sym_graph': True, - 'mg_levels': {'ksp_type': 'gmres', - 'ksp_max_it': 5, - 'pc_type': 'bjacobi', - 'sub_pc_type': 'ilu'}}} - - alpha = 0.51 # off-centering parameter - linear_solver = CompressibleSolver( - eqns, alpha, solver_parameters=params, - overwrite_solver_parameters=True - ) + tau_values = {'rho': 1.0, 'theta': 1.0} + linear_solver = CompressibleSolver(eqns, alpha, tau_values=tau_values) # Time stepper stepper = SemiImplicitQuasiNewton( eqns, io, transported_fields, transport_methods, - linear_solver=linear_solver, alpha=alpha + linear_solver=linear_solver, alpha=alpha, spinup_steps=spinup_steps ) # ------------------------------------------------------------------------ # @@ -189,7 +178,8 @@ def mountain_hydrostatic( bottom_boundary = Constant(exner_surf, domain=mesh) logger.info(f'Solving hydrostatic with bottom Exner of {exner_surf}') compressible_hydrostatic_balance( - eqns, theta_b, rho_b, exner, top=False, exner_boundary=bottom_boundary + eqns, theta_b, rho_b, exner, top=False, exner_boundary=bottom_boundary, + solve_for_rho=True ) # Solve hydrostatic balance again, but now use minimum value from first @@ -243,7 +233,7 @@ def mountain_hydrostatic( theta0.assign(theta_b) rho0.assign(rho_b) - u0.project(as_vector([initial_wind, 0.0]), bcs=eqns.bcs['u']) + u0.project(as_vector([initial_wind, 0.0])) stepper.set_reference_profiles([('rho', rho_b), ('theta', theta_b)]) @@ -268,38 +258,38 @@ def mountain_hydrostatic( '--ncolumns', help="The number of columns in the vertical slice mesh.", type=int, - default=mountain_hydrostatic_defaults['ncolumns'] + default=schaer_mountain_defaults['ncolumns'] ) parser.add_argument( '--nlayers', help="The number of layers for the mesh.", type=int, - default=mountain_hydrostatic_defaults['nlayers'] + default=schaer_mountain_defaults['nlayers'] ) parser.add_argument( '--dt', help="The time step in seconds.", type=float, - default=mountain_hydrostatic_defaults['dt'] + default=schaer_mountain_defaults['dt'] ) parser.add_argument( "--tmax", help="The end time for the simulation in seconds.", type=float, - default=mountain_hydrostatic_defaults['tmax'] + default=schaer_mountain_defaults['tmax'] ) parser.add_argument( '--dumpfreq', help="The frequency at which to dump field output.", type=int, - default=mountain_hydrostatic_defaults['dumpfreq'] + default=schaer_mountain_defaults['dumpfreq'] ) parser.add_argument( '--dirname', help="The name of the directory to write to.", type=str, - default=mountain_hydrostatic_defaults['dirname'] + default=schaer_mountain_defaults['dirname'] ) args, unknown = parser.parse_known_args() - mountain_hydrostatic(**vars(args)) + schaer_mountain(**vars(args)) diff --git a/examples/compressible_euler/skamarock_klemp_hydrostatic.py b/examples/compressible_euler/skamarock_klemp_hydrostatic.py deleted file mode 100644 index f75ad9b10..000000000 --- a/examples/compressible_euler/skamarock_klemp_hydrostatic.py +++ /dev/null @@ -1,204 +0,0 @@ -""" -This example uses the hydrostatic compressible Euler equations to solve the -vertical slice gravity wave test case of Skamarock and Klemp, 1994: -``Efficiency and Accuracy of the Klemp-Wilhelmson Time-Splitting Technique'', -MWR. - -Potential temperature is transported using SUPG, and the degree 1 elements are -used. This also uses a mesh which is one cell thick in the y-direction. -""" - -from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter -from firedrake import ( - as_vector, SpatialCoordinate, PeriodicRectangleMesh, ExtrudedMesh, exp, sin, - Function, pi -) -from gusto import ( - Domain, IO, OutputParameters, SemiImplicitQuasiNewton, SSPRK3, DGUpwind, - TrapeziumRule, SUPGOptions, CourantNumber, Perturbation, - CompressibleParameters, HydrostaticCompressibleEulerEquations, - CompressibleSolver, compressible_hydrostatic_balance -) - -skamarock_klemp_hydrostatic_defaults = { - 'ncolumns': 150, - 'nlayers': 10, - 'dt': 25.0, - 'tmax': 60000., - 'dumpfreq': 1200, - 'dirname': 'skamarock_klemp_hydrostatic' -} - - -def skamarock_klemp_hydrostatic( - ncolumns=skamarock_klemp_hydrostatic_defaults['ncolumns'], - nlayers=skamarock_klemp_hydrostatic_defaults['nlayers'], - dt=skamarock_klemp_hydrostatic_defaults['dt'], - tmax=skamarock_klemp_hydrostatic_defaults['tmax'], - dumpfreq=skamarock_klemp_hydrostatic_defaults['dumpfreq'], - dirname=skamarock_klemp_hydrostatic_defaults['dirname'] -): - - # ------------------------------------------------------------------------ # - # Test case parameters - # ------------------------------------------------------------------------ # - - domain_width = 6.0e6 # Width of domain in x direction (m) - domain_length = 1.0e4 # Length of domain in y direction (m) - domain_height = 1.0e4 # Height of domain (m) - Tsurf = 300. # Temperature at surface (K) - wind_initial = 20. # Initial wind in x direction (m/s) - pert_width = 5.0e3 # Width parameter of perturbation (m) - deltaTheta = 1.0e-2 # Magnitude of theta perturbation (K) - N = 0.01 # Brunt-Vaisala frequency (1/s) - Omega = 0.5e-4 # Planetary rotation rate (1/s) - pressure_gradient_y = -1.0e-4*20 # Prescribed force in y direction (m/s^2) - - # ------------------------------------------------------------------------ # - # Our settings for this set up - # ------------------------------------------------------------------------ # - - element_order = 1 - - # ------------------------------------------------------------------------ # - # Set up model objects - # ------------------------------------------------------------------------ # - - # Domain -- 3D volume mesh - base_mesh = PeriodicRectangleMesh( - ncolumns, 1, domain_width, domain_length, quadrilateral=True - ) - mesh = ExtrudedMesh(base_mesh, nlayers, layer_height=domain_height/nlayers) - domain = Domain(mesh, dt, "RTCF", element_order) - - # Equation - parameters = CompressibleParameters(Omega=Omega) - balanced_pg = as_vector((0., pressure_gradient_y, 0.)) - eqns = HydrostaticCompressibleEulerEquations( - domain, parameters, extra_terms=[("u", balanced_pg)] - ) - - # I/O - output = OutputParameters( - dirname=dirname, dumpfreq=dumpfreq, dump_vtus=True, dump_nc=False, - dumplist=['u'], - ) - diagnostic_fields = [CourantNumber(), Perturbation('theta'), Perturbation('rho')] - io = IO(domain, output, diagnostic_fields=diagnostic_fields) - - # Transport schemes - theta_opts = SUPGOptions() - transported_fields = [ - TrapeziumRule(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, transported_fields, transport_methods, - linear_solver=linear_solver - ) - - # ------------------------------------------------------------------------ # - # Initial conditions - # ------------------------------------------------------------------------ # - - u0 = stepper.fields("u") - rho0 = stepper.fields("rho") - theta0 = stepper.fields("theta") - - # spaces - Vt = domain.spaces("theta") - Vr = domain.spaces("DG") - - # Thermodynamic constants required for setting initial conditions - # and reference profiles - g = parameters.g - - x, _, z = SpatialCoordinate(mesh) - - # N^2 = (g/theta)dtheta/dz => dtheta/dz = theta N^2g => theta=theta_0exp(N^2gz) - thetab = Tsurf*exp(N**2*z/g) - - theta_b = Function(Vt).interpolate(thetab) - rho_b = Function(Vr) - - theta_pert = ( - deltaTheta * sin(pi*z/domain_height) - / (1 + (x - domain_width/2)**2 / pert_width**2) - ) - theta0.interpolate(theta_b + theta_pert) - - compressible_hydrostatic_balance(eqns, theta_b, rho_b, solve_for_rho=True) - - rho0.assign(rho_b) - u0.project(as_vector([wind_initial, 0.0, 0.0])) - - 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=skamarock_klemp_hydrostatic_defaults['ncolumns'] - ) - parser.add_argument( - '--nlayers', - help="The number of layers for the mesh.", - type=int, - default=skamarock_klemp_hydrostatic_defaults['nlayers'] - ) - parser.add_argument( - '--dt', - help="The time step in seconds.", - type=float, - default=skamarock_klemp_hydrostatic_defaults['dt'] - ) - parser.add_argument( - "--tmax", - help="The end time for the simulation in seconds.", - type=float, - default=skamarock_klemp_hydrostatic_defaults['tmax'] - ) - parser.add_argument( - '--dumpfreq', - help="The frequency at which to dump field output.", - type=int, - default=skamarock_klemp_hydrostatic_defaults['dumpfreq'] - ) - parser.add_argument( - '--dirname', - help="The name of the directory to write to.", - type=str, - default=skamarock_klemp_hydrostatic_defaults['dirname'] - ) - args, unknown = parser.parse_known_args() - - skamarock_klemp_hydrostatic(**vars(args)) diff --git a/examples/compressible_euler/skamarock_klemp_nonhydrostatic.py b/examples/compressible_euler/skamarock_klemp_nonhydrostatic.py index 057d34ade..5ec8ab12b 100644 --- a/examples/compressible_euler/skamarock_klemp_nonhydrostatic.py +++ b/examples/compressible_euler/skamarock_klemp_nonhydrostatic.py @@ -4,6 +4,10 @@ ``Efficiency and Accuracy of the Klemp-Wilhelmson Time-Splitting Technique'', MWR. +The domain is smaller than the "hydrostatic" gravity wave test, so that there +is difference between the hydrostatic and non-hydrostatic solutions. The test +can be run with and without a hydrostatic switch. + Potential temperature is transported using SUPG, and the degree 1 elements are used. """ @@ -19,19 +23,20 @@ import numpy as np from gusto import ( Domain, IO, OutputParameters, SemiImplicitQuasiNewton, SSPRK3, DGUpwind, - SUPGOptions, CourantNumber, Perturbation, Gradient, - CompressibleParameters, CompressibleEulerEquations, CompressibleSolver, - compressible_hydrostatic_balance, logger, RichardsonNumber, - RungeKuttaFormulation, SubcyclingOptions + logger, SUPGOptions, Perturbation, CompressibleParameters, + CompressibleEulerEquations, HydrostaticCompressibleEulerEquations, + compressible_hydrostatic_balance, RungeKuttaFormulation, CompressibleSolver, + SubcyclingOptions, hydrostatic_parameters ) skamarock_klemp_nonhydrostatic_defaults = { 'ncolumns': 150, 'nlayers': 10, 'dt': 6.0, - 'tmax': 3600., - 'dumpfreq': 300, - 'dirname': 'skamarock_klemp_nonhydrostatic' + 'tmax': 3000., + 'dumpfreq': 250, + 'dirname': 'skamarock_klemp_nonhydrostatic', + 'hydrostatic': False } @@ -41,7 +46,8 @@ def skamarock_klemp_nonhydrostatic( dt=skamarock_klemp_nonhydrostatic_defaults['dt'], tmax=skamarock_klemp_nonhydrostatic_defaults['tmax'], dumpfreq=skamarock_klemp_nonhydrostatic_defaults['dumpfreq'], - dirname=skamarock_klemp_nonhydrostatic_defaults['dirname'] + dirname=skamarock_klemp_nonhydrostatic_defaults['dirname'], + hydrostatic=skamarock_klemp_nonhydrostatic_defaults['hydrostatic'] ): # ------------------------------------------------------------------------ # @@ -73,18 +79,25 @@ def skamarock_klemp_nonhydrostatic( # Equation parameters = CompressibleParameters() - eqns = CompressibleEulerEquations(domain, parameters) + if hydrostatic: + eqns = HydrostaticCompressibleEulerEquations(domain, parameters) + else: + eqns = CompressibleEulerEquations(domain, parameters) # I/O points_x = np.linspace(0., domain_width, 100) points_z = [domain_height/2.] points = np.array([p for p in itertools.product(points_x, points_z)]) + # Adjust default directory name + if hydrostatic and dirname == skamarock_klemp_nonhydrostatic_defaults['dirname']: + dirname = f'hyd_switch_{dirname}' + # Dumping point data using legacy PointDataOutput is not supported in parallel if COMM_WORLD.size == 1: output = OutputParameters( dirname=dirname, dumpfreq=dumpfreq, pddumpfreq=dumpfreq, - dump_vtus=True, dump_nc=False, + dump_vtus=False, dump_nc=True, point_data=[('theta_perturbation', points)], ) else: @@ -94,14 +107,10 @@ def skamarock_klemp_nonhydrostatic( ) output = OutputParameters( dirname=dirname, dumpfreq=dumpfreq, pddumpfreq=dumpfreq, - dump_vtus=True, dump_nc=True, + dump_vtus=False, dump_nc=True, ) - diagnostic_fields = [ - CourantNumber(), Gradient('u'), Perturbation('theta'), - Gradient('theta_perturbation'), Perturbation('rho'), - RichardsonNumber('theta', parameters.g/Tsurf), Gradient('theta') - ] + diagnostic_fields = [Perturbation('theta')] io = IO(domain, output, diagnostic_fields=diagnostic_fields) # Transport schemes @@ -125,7 +134,13 @@ def skamarock_klemp_nonhydrostatic( ] # Linear solver - linear_solver = CompressibleSolver(eqns) + if hydrostatic: + linear_solver = CompressibleSolver( + eqns, solver_parameters=hydrostatic_parameters, + overwrite_solver_parameters=True + ) + else: + linear_solver = CompressibleSolver(eqns) # Time stepper stepper = SemiImplicitQuasiNewton( @@ -223,6 +238,16 @@ def skamarock_klemp_nonhydrostatic( type=str, default=skamarock_klemp_nonhydrostatic_defaults['dirname'] ) + parser.add_argument( + '--hydrostatic', + help=( + "Whether to use the hydrostatic switch to emulate the " + + "hydrostatic equations. Otherwise use the full non-hydrostatic" + + "equations." + ), + action="store_true", + default=skamarock_klemp_nonhydrostatic_defaults['hydrostatic'] + ) args, unknown = parser.parse_known_args() skamarock_klemp_nonhydrostatic(**vars(args)) diff --git a/examples/compressible_euler/test_compressible_euler_examples.py b/examples/compressible_euler/test_compressible_euler_examples.py index 384824e12..d74d5dade 100644 --- a/examples/compressible_euler/test_compressible_euler_examples.py +++ b/examples/compressible_euler/test_compressible_euler_examples.py @@ -46,66 +46,46 @@ def test_dry_bryan_fritsch_parallel(): test_dry_bryan_fritsch() -# Hydrostatic equations not currently working -@pytest.mark.xfail -def test_mountain_hydrostatic(): - from mountain_hydrostatic import mountain_hydrostatic - test_name = 'mountain_hydrostatic' - mountain_hydrostatic( - ncolumns=20, - nlayers=10, - dt=5.0, - tmax=50.0, - dumpfreq=10, - dirname=make_dirname(test_name) - ) - - -# Hydrostatic equations not currently working -@pytest.mark.xfail -@pytest.mark.parallel(nprocs=4) -def test_mountain_hydrostatic_parallel(): - test_mountain_hydrostatic() - - -# Hydrostatic equations not currently working -@pytest.mark.xfail -def test_skamarock_klemp_hydrostatic(): - from skamarock_klemp_hydrostatic import skamarock_klemp_hydrostatic - test_name = 'skamarock_klemp_hydrostatic' - skamarock_klemp_hydrostatic( +def test_skamarock_klemp_nonhydrostatic(): + from skamarock_klemp_nonhydrostatic import skamarock_klemp_nonhydrostatic + test_name = 'skamarock_klemp_nonhydrostatic' + skamarock_klemp_nonhydrostatic( ncolumns=30, nlayers=5, dt=6.0, tmax=60.0, dumpfreq=10, - dirname=make_dirname(test_name) + dirname=make_dirname(test_name), + hydrostatic=False ) -# Hydrostatic equations not currently working -@pytest.mark.xfail @pytest.mark.parallel(nprocs=2) -def test_skamarock_klemp_hydrostatic_parallel(): - test_skamarock_klemp_hydrostatic() +def test_skamarock_klemp_nonhydrostatic_parallel(): + test_skamarock_klemp_nonhydrostatic() -def test_skamarock_klemp_nonhydrostatic(): +# Hydrostatic equations not currently working +@pytest.mark.xfail +def test_hyd_switch_skamarock_klemp_nonhydrostatic(): from skamarock_klemp_nonhydrostatic import skamarock_klemp_nonhydrostatic - test_name = 'skamarock_klemp_nonhydrostatic' + test_name = 'hyd_switch_skamarock_klemp_nonhydrostatic' skamarock_klemp_nonhydrostatic( ncolumns=30, nlayers=5, dt=6.0, tmax=60.0, dumpfreq=10, - dirname=make_dirname(test_name) + dirname=make_dirname(test_name), + hydrostatic=True ) +# Hydrostatic equations not currently working +@pytest.mark.xfail @pytest.mark.parallel(nprocs=2) -def test_skamarock_klemp_nonhydrostatic_parallel(): - test_skamarock_klemp_nonhydrostatic() +def test_hyd_switch_skamarock_klemp_nonhydrostatic_parallel(): + test_hyd_switch_skamarock_klemp_nonhydrostatic() def test_straka_bubble(): diff --git a/figures/compressible_euler/schaer_mountain_final.png b/figures/compressible_euler/schaer_mountain_final.png new file mode 100644 index 000000000..a5289cca3 Binary files /dev/null and b/figures/compressible_euler/schaer_mountain_final.png differ diff --git a/figures/compressible_euler/schaer_mountain_initial.png b/figures/compressible_euler/schaer_mountain_initial.png new file mode 100644 index 000000000..5e78c469a Binary files /dev/null and b/figures/compressible_euler/schaer_mountain_initial.png differ diff --git a/figures/compressible_euler/skamarock_klemp_nonhydrostatic_final.png b/figures/compressible_euler/skamarock_klemp_nonhydrostatic_final.png index c1187a1cd..e76a82582 100644 Binary files a/figures/compressible_euler/skamarock_klemp_nonhydrostatic_final.png and b/figures/compressible_euler/skamarock_klemp_nonhydrostatic_final.png differ diff --git a/gusto/solvers/linear_solvers.py b/gusto/solvers/linear_solvers.py index ca684c6af..213d09df2 100644 --- a/gusto/solvers/linear_solvers.py +++ b/gusto/solvers/linear_solvers.py @@ -55,7 +55,7 @@ def __init__(self, equations, alpha=0.5, tau_values=None, """ self.equations = equations self.dt = equations.domain.dt - self.alpha = alpha + self.alpha = Constant(alpha) self.tau_values = tau_values if tau_values is not None else {} if solver_parameters is not None: @@ -811,9 +811,10 @@ def __init__(self, equation, alpha, reference_dependent=True): lambda t: Term(t.get(linearisation).form, t.labels), drop) + self.alpha = Constant(alpha) dt = equation.domain.dt W = equation.function_space - beta = dt*alpha + beta = dt*self.alpha # Split up the rhs vector (symbolically) self.xrhs = Function(W) diff --git a/gusto/solvers/parameters.py b/gusto/solvers/parameters.py index a7f61edb0..f7f0d93dd 100644 --- a/gusto/solvers/parameters.py +++ b/gusto/solvers/parameters.py @@ -4,7 +4,7 @@ """ from gusto.core.function_spaces import is_cg -__all__ = ['mass_parameters'] +__all__ = ['mass_parameters', 'hydrostatic_parameters'] def mass_parameters(V, spaces=None, ignore_vertical=True): @@ -97,3 +97,28 @@ def mass_parameters(V, spaces=None, ignore_vertical=True): } return parameters + + +hydrostatic_parameters = { + 'mat_type': 'matfree', + 'ksp_type': 'preonly', + 'pc_type': 'python', + 'pc_python_type': 'firedrake.SCPC', + # Velocity mass operator is singular in the hydrostatic case. + # So for reconstruction, we eliminate rho into u + 'pc_sc_eliminate_fields': '1, 0', + 'condensed_field': { + 'ksp_type': 'fgmres', + 'ksp_rtol': 1.0e-8, + 'ksp_atol': 1.0e-8, + 'ksp_max_it': 100, + 'pc_type': 'gamg', + 'pc_gamg_sym_graph': True, + 'mg_levels': { + 'ksp_type': 'gmres', + 'ksp_max_it': 5, + 'pc_type': 'bjacobi', + 'sub_pc_type': 'ilu' + } + } +} diff --git a/gusto/time_discretisation/time_discretisation.py b/gusto/time_discretisation/time_discretisation.py index 3e50be6b7..db71412ca 100644 --- a/gusto/time_discretisation/time_discretisation.py +++ b/gusto/time_discretisation/time_discretisation.py @@ -418,11 +418,12 @@ def update_subcycling(self): # Cap number of subcycles if max_subcycles is not None: - self.ncycles = min(self.ncycles, max_subcycles) - logger.warning( - 'Adaptive subcycling: capping number of subcycles at ' - f'{max_subcycles}' - ) + if self.ncycles > max_subcycles: + logger.warning( + 'Adaptive subcycling: capping number of subcycles at ' + f'{max_subcycles}' + ) + self.ncycles = max_subcycles logger.debug(f'Performing {self.ncycles} subcycles') self.dt.assign(self.original_dt/self.ncycles) diff --git a/gusto/timestepping/semi_implicit_quasi_newton.py b/gusto/timestepping/semi_implicit_quasi_newton.py index 19142dbaa..deacb6af2 100644 --- a/gusto/timestepping/semi_implicit_quasi_newton.py +++ b/gusto/timestepping/semi_implicit_quasi_newton.py @@ -39,7 +39,8 @@ def __init__(self, equation_set, io, transport_schemes, spatial_methods, slow_physics_schemes=None, fast_physics_schemes=None, alpha=Constant(0.5), off_centred_u=False, num_outer=2, num_inner=2, accelerator=False, - predictor=None, reference_update_freq=None): + predictor=None, reference_update_freq=None, + spinup_steps=0): """ Args: equation_set (:class:`PrognosticEquationSet`): the prognostic @@ -106,15 +107,24 @@ def __init__(self, equation_set, io, transport_schemes, spatial_methods, time step. Setting it to None turns off the update, and reference profiles will remain at their initial values. Defaults to None. + spinup_steps (int, optional): the number of steps to run the model + in "spin-up" mode, where the alpha parameter is set to 1.0. + Defaults to 0, which corresponds to no spin-up. """ self.num_outer = num_outer self.num_inner = num_inner - self.alpha = alpha + self.alpha = Constant(alpha) self.predictor = predictor self.accelerator = accelerator + + # Options relating to reference profiles and spin-up + self._alpha_original = Constant(alpha) self.reference_update_freq = reference_update_freq self.to_update_ref_profile = False + self.spinup_steps = spinup_steps + self.spinup_begun = False + self.spinup_done = False # Flag for if we have simultaneous transport self.simult = False @@ -353,13 +363,43 @@ def update_reference_profiles(self): if float(self.t) + self.reference_update_freq > self.last_ref_update_time: self.equation.X_ref.assign(self.x.n(self.field_name)) self.last_ref_update_time = float(self.t) - if hasattr(self.linear_solver, 'update_reference_profiles'): - self.linear_solver.update_reference_profiles() + self.linear_solver.update_reference_profiles() elif self.to_update_ref_profile: - if hasattr(self.linear_solver, 'update_reference_profiles'): - self.linear_solver.update_reference_profiles() - self.to_update_ref_profile = False + self.linear_solver.update_reference_profiles() + self.to_update_ref_profile = False + + def start_spinup(self): + """ + Initialises the spin-up period, so that the scheme is implicit by + setting the off-centering parameter alpha to be 1. + """ + logger.debug('Starting spin-up period') + # Update alpha + self.alpha.assign(1.0) + self.linear_solver.alpha.assign(1.0) + # We need to tell solvers that they may need rebuilding + self.linear_solver.update_reference_profiles() + self.forcing.solvers['explicit'].invalidate_jacobian() + self.forcing.solvers['implicit'].invalidate_jacobian() + # This only needs doing once, so update the flag + self.spinup_begun = True + + def finish_spinup(self): + """ + Finishes the spin-up period, returning the off-centering parameter + to its original value. + """ + logger.debug('Finishing spin-up period') + # Update alpha + self.alpha.assign(self._alpha_original) + self.linear_solver.alpha.assign(self._alpha_original) + # We need to tell solvers that they may need rebuilding + self.linear_solver.update_reference_profiles() + self.forcing.solvers['explicit'].invalidate_jacobian() + self.forcing.solvers['implicit'].invalidate_jacobian() + # This only needs doing once, so update the flag + self.spinup_done = True def timestep(self): """Defines the timestep""" @@ -376,6 +416,13 @@ def timestep(self): # Update reference profiles -------------------------------------------- self.update_reference_profiles() + # Are we in spin-up period? -------------------------------------------- + # Note: steps numbered from 1 onwards + if self.step < self.spinup_steps + 1 and not self.spinup_begun: + self.start_spinup() + elif self.step >= self.spinup_steps + 1 and not self.spinup_done: + self.finish_spinup() + # Slow physics --------------------------------------------------------- x_after_slow(self.field_name).assign(xn(self.field_name)) if len(self.slow_physics_schemes) > 0: @@ -526,7 +573,7 @@ def __init__(self, equation, alpha): map_if_false=drop) # the explicit forms are multiplied by (1-alpha) and moved to the rhs - L_explicit = -(1-alpha)*dt*residual.label_map( + L_explicit = -(Constant(1)-alpha)*dt*residual.label_map( lambda t: any(t.has_label(time_derivative, hydrostatic, *implicit_terms, return_tuple=True)), diff --git a/plotting/compressible_euler/plot_schaer_mountain.py b/plotting/compressible_euler/plot_schaer_mountain.py new file mode 100644 index 000000000..f3d3f7e3d --- /dev/null +++ b/plotting/compressible_euler/plot_schaer_mountain.py @@ -0,0 +1,181 @@ +""" +Plots the Schär mountain test case. + +This plots: +(a) w @ t = 5 hr, (b) theta perturbation @ t = 5 hr +""" +from os.path import abspath, dirname +import matplotlib.pyplot as plt +from netCDF4 import Dataset +import numpy as np +from tomplot import ( + set_tomplot_style, tomplot_cmap, plot_contoured_field, + add_colorbar_ax, tomplot_field_title, tomplot_contours, + extract_gusto_coords, extract_gusto_field, reshape_gusto_data +) + +test = 'schaer_mountain' + +# ---------------------------------------------------------------------------- # +# Directory for results and plots +# ---------------------------------------------------------------------------- # +# When copying this example these paths need editing, which will usually involve +# removing the abspath part to set directory paths relative to this file +results_file_name = f'{abspath(dirname(__file__))}/../../results/{test}/field_output.nc' +plot_stem = f'{abspath(dirname(__file__))}/../../figures/compressible_euler/{test}' + +# ---------------------------------------------------------------------------- # +# Final plot details +# ---------------------------------------------------------------------------- # +final_field_names = ['u_z', 'theta_perturbation', 'u_z', 'theta_perturbation'] +final_colour_schemes = ['PiYG', 'RdBu_r', 'PiYG', 'RdBu_r'] +final_field_labels = [ + r'$w$ (m s$^{-1}$)', r'$\Delta\theta$ (K)', + r'$w$ (m s$^{-1}$)', r'$\Delta\theta$ (K)' +] +final_contours = [ + np.linspace(-1.1, 1.1, 23), np.linspace(-1.4, 1.4, 15), + np.linspace(-1.1, 1.1, 23), np.linspace(-1.4, 1.4, 31) +] + +# ---------------------------------------------------------------------------- # +# Initial plot details +# ---------------------------------------------------------------------------- # +initial_field_names = ['Exner', 'theta'] +initial_colour_schemes = ['PuBu', 'Reds'] +initial_field_labels = [r'$\Pi$', r'$\theta$ (K)'] + +# ---------------------------------------------------------------------------- # +# General options +# ---------------------------------------------------------------------------- # +contour_method = 'contour' # Need to use this method to show mountains! +xlims = [0., 100.] +ylims = [0., 30.] + +# Things that are likely the same for all plots -------------------------------- +set_tomplot_style() +data_file = Dataset(results_file_name, 'r') + +# ---------------------------------------------------------------------------- # +# INITIAL PLOTTING +# ---------------------------------------------------------------------------- # + +fig, axarray = plt.subplots(1, 2, figsize=(18, 6), sharex='all', sharey='all') +time_idx = 0 + +for i, (ax, field_name, colour_scheme, field_label) in \ + enumerate(zip( + axarray.flatten(), initial_field_names, initial_colour_schemes, + initial_field_labels + )): + + # Data extraction ---------------------------------------------------------- + field_data = extract_gusto_field(data_file, field_name, time_idx=time_idx) + coords_X, coords_Y = extract_gusto_coords(data_file, field_name) + field_data, coords_X, coords_Y = \ + reshape_gusto_data(field_data, coords_X, coords_Y) + time = data_file['time'][time_idx] + + contours = tomplot_contours(field_data) + cmap, lines = tomplot_cmap(contours, colour_scheme) + + # Plot data ---------------------------------------------------------------- + cf, _ = plot_contoured_field( + ax, coords_X, coords_Y, field_data, contour_method, contours, + cmap=cmap, line_contours=lines + ) + + add_colorbar_ax(ax, cf, field_label, location='bottom') + tomplot_field_title(ax, None, minmax=True, field_data=field_data) + + # Labels ------------------------------------------------------------------- + if i == 0: + ax.set_ylabel(r'$z$ (km)', labelpad=-20) + ax.set_ylim(ylims) + ax.set_yticks(ylims) + ax.set_yticklabels(ylims) + + ax.set_xlabel(r'$x$ (km)', labelpad=-10) + ax.set_xlim(xlims) + ax.set_xticks(xlims) + ax.set_xticklabels(xlims) + +# Save figure ------------------------------------------------------------------ +fig.suptitle(f't = {time:.1f} s') +fig.subplots_adjust(wspace=0.25) +plot_name = f'{plot_stem}_initial.png' +print(f'Saving figure to {plot_name}') +fig.savefig(plot_name, bbox_inches='tight') +plt.close() + +# ---------------------------------------------------------------------------- # +# FINAL PLOTTING +# ---------------------------------------------------------------------------- # +xlims_zoom = [30., 70.] +ylims_zoom = [0., 12.] + +fig, axarray = plt.subplots(2, 2, figsize=(18, 12), sharex='row', sharey='row') +time_idx = 1 + +for i, (ax, field_name, colour_scheme, field_label, contours) in \ + enumerate(zip( + axarray.flatten(), final_field_names, final_colour_schemes, + final_field_labels, final_contours + )): + + # Data extraction ---------------------------------------------------------- + field_data = extract_gusto_field(data_file, field_name, time_idx=time_idx) + coords_X, coords_Y = extract_gusto_coords(data_file, field_name) + time = data_file['time'][time_idx] + + field_data, coords_X, coords_Y = \ + reshape_gusto_data(field_data, coords_X, coords_Y) + + cmap, lines = tomplot_cmap(contours, colour_scheme, remove_contour=0.0) + + # Plot data ---------------------------------------------------------------- + cf, _ = plot_contoured_field( + ax, coords_X, coords_Y, field_data, contour_method, contours, + cmap=cmap, line_contours=lines + ) + + add_colorbar_ax(ax, cf, field_label, location='bottom') + + if i in [0, 1]: + # Only print min/max for top plots + tomplot_field_title( + ax, None, minmax=True, field_data=field_data, minmax_format='.3f' + ) + + # Labels ------------------------------------------------------------------- + ax.set_xlabel(r'$x$ (km)', labelpad=-10) + + if i in [0, 1]: + ax.set_xlim(xlims) + ax.set_ylim(ylims) + ax.set_xticks(xlims) + ax.set_xticklabels(xlims) + else: + ax.set_xlim(xlims_zoom) + ax.set_ylim(ylims_zoom) + ax.set_xticks(xlims_zoom) + ax.set_xticklabels(xlims_zoom) + + if i == 0: + ax.set_ylabel(r'$z$ (km)', labelpad=-20) + ax.set_yticks(ylims) + ax.set_yticklabels(ylims) + + elif i == 2: + ax.set_ylabel(r'$z$ (km)', labelpad=-20) + ax.set_yticks(ylims_zoom) + ax.set_yticklabels(ylims_zoom) + + +# Save figure ------------------------------------------------------------------ +fig.suptitle(f't = {time:.1f} s') +fig.subplots_adjust(wspace=0.25) +plot_name = f'{plot_stem}_final.png' +print(f'Saving figure to {plot_name}') +fig.savefig(plot_name, bbox_inches='tight') +plt.close() diff --git a/plotting/compressible_euler/plot_skamarock_klemp_nonhydrostatic.py b/plotting/compressible_euler/plot_skamarock_klemp_nonhydrostatic.py index 668d75dc8..566948b41 100644 --- a/plotting/compressible_euler/plot_skamarock_klemp_nonhydrostatic.py +++ b/plotting/compressible_euler/plot_skamarock_klemp_nonhydrostatic.py @@ -3,7 +3,7 @@ This plots the initial conditions @ t = 0 s, with (a) theta perturbation, (b) theta -and the final state @ t = 3600 s, with +and the final state @ t = 3000 s, with (a) theta perturbation, (b) a 1D slice through the wave """ @@ -11,6 +11,7 @@ import matplotlib.pyplot as plt import numpy as np from netCDF4 import Dataset +import pandas as pd from tomplot import ( set_tomplot_style, tomplot_cmap, plot_contoured_field, add_colorbar_ax, tomplot_field_title, extract_gusto_coords, @@ -49,7 +50,6 @@ # General options # ---------------------------------------------------------------------------- # contour_method = 'tricontour' -xlims = [0, 300.0] ylims = [0, 10.0] # Things that are likely the same for all plots -------------------------------- @@ -59,6 +59,8 @@ # ---------------------------------------------------------------------------- # # INITIAL PLOTTING # ---------------------------------------------------------------------------- # +xlims = [0, 300.0] + fig, axarray = plt.subplots(1, 2, figsize=(12, 6), sharex='all', sharey='all') time_idx = 0 @@ -107,6 +109,9 @@ # ---------------------------------------------------------------------------- # # FINAL PLOTTING # ---------------------------------------------------------------------------- # +x_offset = -3000.0*20/1000.0 +xlims = [-x_offset, 300.0-x_offset] + fig, axarray = plt.subplots(2, 1, figsize=(8, 8), sharex='all') time_idx = -1 @@ -115,6 +120,21 @@ coords_X, coords_Y = extract_gusto_coords(data_file, final_field_name) time = data_file['time'][time_idx] +# Wave has wrapped around periodic boundary, so shift the coordinates +coords_X = np.where(coords_X < xlims[0], coords_X + 300.0, coords_X) + +# Sort data given the change in coordinates +data_dict = { + 'X': coords_X, + 'Y': coords_Y, + 'field': field_data +} +data_frame = pd.DataFrame(data_dict) +data_frame.sort_values(by=['X', 'Y'], inplace=True) +coords_X = data_frame['X'].values[:] +coords_Y = data_frame['Y'].values[:] +field_data = data_frame['field'].values[:] + # Plot 2D data ----------------------------------------------------------------- ax = axarray[0]