New model

examples/compressible_euler/dry_bryan_fritsch.py

@@ -13,8 +13,8 @@
     LinearVariationalProblem, LinearVariationalSolver
 )
 from gusto import (
-    Domain, IO, OutputParameters, SemiImplicitQuasiNewton, SSPRK3, DGUpwind,
-    RecoverySpaces, BoundaryMethod, Perturbation, CompressibleParameters,
+    OutputParameters, LowestOrderModel,
+    Perturbation, CompressibleParameters,
     CompressibleEulerEquations, compressible_hydrostatic_balance
 )
 
@@ -52,7 +52,6 @@ def dry_bryan_fritsch(
     # Our settings for this set up
     # ------------------------------------------------------------------------ #
 
-    element_order = 0
     u_eqn_type = 'vector_advection_form'
 
     # ------------------------------------------------------------------------ #
@@ -62,61 +61,37 @@ def dry_bryan_fritsch(
     # Domain
     base_mesh = IntervalMesh(ncolumns, domain_width)
     mesh = ExtrudedMesh(base_mesh, nlayers, layer_height=domain_height/nlayers)
-    domain = Domain(mesh, dt, "CG", element_order)
 
     # Equation
     params = CompressibleParameters(mesh)
-    eqns = CompressibleEulerEquations(
-        domain, params, u_transport_option=u_eqn_type,
-        no_normal_flow_bc_ids=[1, 2]
-    )
+    eqns = CompressibleEulerEquations
+
+    # Model
+    model = LowestOrderModel(mesh, dt, params, eqns,
+                             u_transport_option=u_eqn_type,
+                             family="CG", no_normal_flow_bc_ids=[1, 2])
 
     # I/O
     output = OutputParameters(
         dirname=dirname, dumpfreq=dumpfreq, dump_vtus=False, dump_nc=True,
         dumplist=['rho']
     )
     diagnostic_fields = [Perturbation('theta')]
-    io = IO(domain, output, diagnostic_fields=diagnostic_fields)
-
-    # Transport schemes -- set up options for using recovery wrapper
-    boundary_methods = {'DG': BoundaryMethod.taylor}
-
-    recovery_spaces = RecoverySpaces(
-        domain, boundary_methods, use_vector_spaces=True
-    )
 
-    u_opts = recovery_spaces.HDiv_options
-    rho_opts = recovery_spaces.DG_options
-    theta_opts = recovery_spaces.theta_options
-
-    transported_fields = [
-        SSPRK3(domain, "rho", options=rho_opts),
-        SSPRK3(domain, "theta", options=theta_opts),
-        SSPRK3(domain, "u", options=u_opts)
-    ]
-
-    transport_methods = [
-        DGUpwind(eqns, field) for field in ["u", "rho", "theta"]
-    ]
-
-    # Time stepper
-    stepper = SemiImplicitQuasiNewton(
-        eqns, io, transported_fields, transport_methods,
-        tau_values={'rho': 1.0, 'theta': 1.0}
-    )
+    model.setup(output, diagnostic_fields=diagnostic_fields)
 
     # ------------------------------------------------------------------------ #
     # Initial conditions
     # ------------------------------------------------------------------------ #
 
+    stepper = model.stepper
     u0 = stepper.fields("u")
     rho0 = stepper.fields("rho")
     theta0 = stepper.fields("theta")
 
     # spaces
-    Vt = domain.spaces("theta")
-    Vr = domain.spaces("DG")
+    Vt = model.domain.spaces("theta")
+    Vr = model.domain.spaces("DG")
     x, z = SpatialCoordinate(mesh)
 
     # Define constant theta_e and water_t
@@ -127,7 +102,7 @@ def dry_bryan_fritsch(
     u0.project(as_vector([zero, zero]))
 
     # Calculate hydrostatic fields
-    compressible_hydrostatic_balance(eqns, theta_b, rho0, solve_for_rho=True)
+    compressible_hydrostatic_balance(model.equation, theta_b, rho0, solve_for_rho=True)
 
     # make mean fields
     rho_b = Function(Vr).assign(rho0)
@@ -161,7 +136,7 @@ def dry_bryan_fritsch(
     # Run
     # ------------------------------------------------------------------------ #
 
-    stepper.run(t=0, tmax=tmax)
+    model.run(t=0, tmax=tmax)
 
 # ---------------------------------------------------------------------------- #
 # MAIN

examples/compressible_euler/schaer_mountain.py

@@ -15,11 +15,11 @@
     SpatialCoordinate, exp, pi, cos, Function, Mesh, Constant
 )
 from gusto import (
-    Domain, CompressibleParameters, logger,
-    OutputParameters, IO, SSPRK3, DGUpwind, SemiImplicitQuasiNewton,
+    CompressibleParameters, logger,
+    OutputParameters, Model,
     compressible_hydrostatic_balance, SpongeLayerParameters, Exner, ZComponent,
-    Perturbation, SUPGOptions, TrapeziumRule, MaxKernel, MinKernel,
-    CompressibleEulerEquations, SubcyclingOptions, RungeKuttaFormulation
+    Perturbation, MaxKernel, MinKernel,
+    CompressibleEulerEquations, SubcyclingOptions
 )
 
 schaer_mountain_defaults = {
@@ -66,7 +66,6 @@ def schaer_mountain(
 
     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'
 
     # ------------------------------------------------------------------------ #
@@ -93,16 +92,19 @@ def schaer_mountain(
     new_coords = Function(Vc).interpolate(xexpr)
     mesh = Mesh(new_coords)
     mesh._base_mesh = base_mesh  # Force new mesh to inherit original base mesh
-    domain = Domain(mesh, dt, "CG", element_order)
 
     # Equation
     parameters = CompressibleParameters(mesh, g=g, cp=cp)
     sponge_params = SpongeLayerParameters(
         mesh, H=domain_height, z_level=domain_height-sponge_depth, mubar=mu_dt/dt
     )
-    eqns = CompressibleEulerEquations(
-        domain, parameters, sponge_options=sponge_params, u_transport_option=u_eqn_type
-    )
+    eqns = CompressibleEulerEquations
+
+    # Model
+    model = Model(mesh, dt, parameters, eqns,
+                  u_transport_option=u_eqn_type,
+                  sponge_options=sponge_params,
+                  family="CG")
 
     # I/O
     output = OutputParameters(
@@ -112,46 +114,25 @@ def schaer_mountain(
         Exner(parameters), ZComponent('u'), Perturbation('theta'),
         Perturbation('rho')
     ]
-    io = IO(domain, output, diagnostic_fields=diagnostic_fields)
 
-    # Transport schemes
+    # Setup model
     subcycling_opts = SubcyclingOptions(subcycle_by_courant=0.25)
-    theta_opts = SUPGOptions()
-    transported_fields = [
-        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", advective_then_flux=True),
-        DGUpwind(eqns, "theta", ibp=theta_opts.ibp)
-    ]
-
-    # Time stepper
-    tau_values = {'rho': 1.0, 'theta': 1.0}
-    stepper = SemiImplicitQuasiNewton(
-        eqns, io, transported_fields, transport_methods,
-        alpha=alpha, tau_values=tau_values, spinup_steps=spinup_steps
-    )
+    model.setup(output, subcycling_options=subcycling_opts,
+                diagnostic_fields=diagnostic_fields,
+                alpha=alpha, spinup_steps=spinup_steps)
 
     # ------------------------------------------------------------------------ #
     # Initial conditions
     # ------------------------------------------------------------------------ #
 
+    stepper = model.stepper
     u0 = stepper.fields("u")
     rho0 = stepper.fields("rho")
     theta0 = stepper.fields("theta")
 
     # spaces
-    Vt = domain.spaces("theta")
-    Vr = domain.spaces("DG")
+    Vt = model.domain.spaces("theta")
+    Vr = model.domain.spaces("DG")
 
     # Thermodynamic constants required for setting initial conditions
     # and reference profiles
@@ -175,8 +156,8 @@ def schaer_mountain(
     bottom_boundary = Constant(exner_surf)
     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,
-        solve_for_rho=True
+        model.equation, 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
@@ -185,7 +166,8 @@ def schaer_mountain(
     top_boundary = Constant(top_value)
     logger.info(f'Solving hydrostatic with top Exner of {top_value}')
     compressible_hydrostatic_balance(
-        eqns, theta_b, rho_b, exner, top=True, exner_boundary=top_boundary
+        model.equation, theta_b, rho_b, exner, top=True,
+        exner_boundary=top_boundary
     )
 
     max_bottom_value = max_kernel.apply(exner)
@@ -209,7 +191,8 @@ def schaer_mountain(
         )
 
         compressible_hydrostatic_balance(
-            eqns, theta_b, rho_b, exner, top=True, exner_boundary=top_boundary
+            model.equation, theta_b, rho_b, exner, top=True,
+            exner_boundary=top_boundary
         )
 
         max_bottom_value = max_kernel.apply(exner)
@@ -224,8 +207,8 @@ def schaer_mountain(
 
     # Perform a final solve to obtain hydrostatically balanced rho
     compressible_hydrostatic_balance(
-        eqns, theta_b, rho_b, exner, top=True, exner_boundary=top_boundary,
-        solve_for_rho=True
+        model.equation, theta_b, rho_b, exner, top=True,
+        exner_boundary=top_boundary, solve_for_rho=True
     )
 
     theta0.assign(theta_b)
@@ -238,7 +221,7 @@ def schaer_mountain(
     # Run
     # ------------------------------------------------------------------------ #
 
-    stepper.run(t=0, tmax=tmax)
+    model.run(t=0, tmax=tmax)
 
 # ---------------------------------------------------------------------------- #
 # MAIN