From 89a7d9b4b05999fdbdf0ad97cb9c307b503f5466 Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Tue, 23 Jul 2024 18:19:51 +0100 Subject: [PATCH 1/4] interface to create form_mass/function from gusto Equation --- case_studies/gusto/advection.py | 171 ++++++++++++++++++++++++++ case_studies/gusto/diffusion.py | 66 ++++++++++ case_studies/gusto/gusto_interface.py | 66 ++++++++++ case_studies/gusto/lswe.py | 92 ++++++++++++++ 4 files changed, 395 insertions(+) create mode 100644 case_studies/gusto/advection.py create mode 100644 case_studies/gusto/diffusion.py create mode 100644 case_studies/gusto/gusto_interface.py create mode 100644 case_studies/gusto/lswe.py diff --git a/case_studies/gusto/advection.py b/case_studies/gusto/advection.py new file mode 100644 index 00000000..325f2f39 --- /dev/null +++ b/case_studies/gusto/advection.py @@ -0,0 +1,171 @@ + +from firedrake import (SpatialCoordinate, PeriodicIntervalMesh, exp, as_vector, + norm, Constant, conditional, sqrt, VectorFunctionSpace) +from firedrake.petsc import PETSc +from gusto import * + +from utils.serial import SerialMiniApp + +Print = PETSc.Sys.Print + + +def asq_forms(equation): + from firedrake.fml import ( + Term, all_terms, drop, + replace_subject, replace_test_function + ) + from gusto.labels import time_derivative, prognostic + NullTerm = Term(None) + + V = equation.function_space + ncpts = len(V) + residual = equation.residual + + # split equation into mass matrix and linear operator + mass = residual.label_map( + lambda t: t.has_label(time_derivative), + map_if_false=drop) + + function = residual.label_map( + lambda t: t.has_label(time_derivative), + map_if_true=drop) + + # generate ufl for mass matrix over given trial/tests + def form_mass(*trials_and_tests): + trials = trials_and_tests[:ncpts] + tests = trials_and_tests[ncpts:] + m = mass.label_map( + all_terms, + replace_test_function(tests[0])) + m = m.label_map( + all_terms, + replace_subject(trials[0])) + return m.form + + # generate ufl for linear operator over given trial/tests + def form_function(*trials_and_tests): + trials = trials_and_tests[:ncpts] + tests = trials_and_tests[ncpts:] + + f = NullTerm + for i in range(ncpts): + fi = function.label_map( + lambda t: t.get(prognostic) == equation.field_name, + lambda t: Term( + split_form(t.form)[i].form, + t.labels), + map_if_false=drop) + + fi = fi.label_map( + all_terms, + replace_test_function(tests[i])) + + fi = fi.label_map( + all_terms, + replace_subject(trials[i])) + + f += fi + f = f.label_map(lambda t: t is NullTerm, drop) + return f.form + + return form_mass, form_function + + +# Set up model objects +# ------------------------------------------------------------------------ # + +# Domain +dt = 0.02 +tmax = 1.0 +L = 10 +mesh = PeriodicIntervalMesh(20, L) +domain = Domain(mesh, dt, "CG", 1) + +# Equation +diffusion_params = DiffusionParameters(kappa=0.75, mu=5) +V = domain.spaces("DG") +Vu = VectorFunctionSpace(mesh, "CG", 1) + +equation = AdvectionDiffusionEquation(domain, V, "f", Vu=Vu, + diffusion_parameters=diffusion_params) +spatial_methods = [DGUpwind(equation, "f"), + InteriorPenaltyDiffusion(equation, "f", diffusion_params)] + +# I/O +output = OutputParameters(dirname=str('tmp'), dumpfreq=25) +io = IO(domain, output) + +# Time stepper +stepper = PrescribedTransport(equation, TrapeziumRule(domain), io, spatial_methods) + +# ------------------------------------------------------------------------ # +# Initial conditions +# ------------------------------------------------------------------------ # + +x = SpatialCoordinate(mesh) +xc_init = 0.25*L +xc_end = 0.75*L +umax = 0.5*L/tmax +# umax = 1 + +# Get minimum distance on periodic interval to xc +x_init = conditional(sqrt((x[0] - xc_init)**2) < 0.5*L, + x[0] - xc_init, L + x[0] - xc_init) + +x_end = conditional(sqrt((x[0] - xc_end)**2) < 0.5*L, + x[0] - xc_end, L + x[0] - xc_end) + +f_init = 5.0 +f_end = f_init / 2.0 +f_width_init = L / 10.0 +f_width_end = f_width_init * 2.0 +f_init_expr = f_init*exp(-(x_init / f_width_init)**2) +f_end_expr = f_end*exp(-(x_end / f_width_end)**2) + +stepper.fields('f').interpolate(f_init_expr) +stepper.fields('u').interpolate(as_vector([Constant(umax)])) +f_end = stepper.fields('f_end', space=V) +f_end.interpolate(f_end_expr) + +# ------------------------------------------------------------------------ # +# Run Gusto stepper +# ------------------------------------------------------------------------ # + +stepper.run(0, tmax=tmax) + +# ------------------------------------------------------------------------ # +# Set up asQ stepper +# ------------------------------------------------------------------------ # + +Print(f"{equation.function_space = }") +Print(f"{V = }") +params = { + 'snes_type': 'ksponly', + 'ksp_type': 'preonly', + 'pc_type': 'lu', + 'ksp_converged_reason': None, +} +form_mass, form_function = asq_forms(equation) +theta = 0.5 +f0 = Function(V).interpolate(f_init_expr) +asqstep = SerialMiniApp(dt, theta, f0, form_mass, form_function, + solver_parameters=params) + +# ------------------------------------------------------------------------ # +# Run Gusto stepper +# ------------------------------------------------------------------------ # + +asqstep.solve(nt=stepper.step-1) + +# ------------------------------------------------------------------------ # +# Check errors +# ------------------------------------------------------------------------ # + +gusto_error = norm(stepper.fields('f') - f_end) / norm(f_end) +Print(f"{gusto_error = }") + +asq_error = norm(asqstep.w1 - f_end) / norm(f_end) +Print(f"{asq_error = }") + +gusto_asq_error = norm(stepper.fields('f') - asqstep.w1) / norm(f_end) +Print(f"{gusto_asq_error = }") diff --git a/case_studies/gusto/diffusion.py b/case_studies/gusto/diffusion.py new file mode 100644 index 00000000..5d55d15e --- /dev/null +++ b/case_studies/gusto/diffusion.py @@ -0,0 +1,66 @@ +from firedrake import ( + PeriodicSquareMesh, SpatialCoordinate, + exp, errornorm, Function) +from gusto import * +from firedrake.petsc import PETSc +Print = PETSc.Sys.Print + +from utils.serial import SerialMiniApp +from gusto_interface import asq_forms + +dt = 0.01 +L = 10. +mesh = PeriodicSquareMesh(10, 10, L, direction='x') + +output = OutputParameters(dirname=str('output'), dumpfreq=25) +domain = Domain(mesh, dt, family="CG", degree=1) +io = IO(domain, output) + +tmax = 1.0 +x = SpatialCoordinate(mesh) +f_init = exp(-((x[0]-0.5*L)**2 + (x[1]-0.5*L)**2)) + +tol = 5.e-2 +kappa = 1. + +f_end_expr = (1/(1+4*tmax))*f_init**(1/(1+4*tmax)) + +V = domain.spaces("DG") +Print(f"{V = }") + +mu = 5. + +# gusto stepper + +diffusion_params = DiffusionParameters(kappa=kappa, mu=mu) +eqn = DiffusionEquation(domain, V, "f", diffusion_parameters=diffusion_params) +diffusion_scheme = BackwardEuler(domain) +diffusion_methods = [InteriorPenaltyDiffusion(eqn, "f", diffusion_params)] +timestepper = Timestepper(eqn, diffusion_scheme, io, spatial_methods=diffusion_methods) +timestepper.fields("f").interpolate(f_init) + +# asq stepper +form_mass, form_function = asq_forms(eqn) +theta = 1 +f0 = Function(V).interpolate(f_init) +params = { + 'snes_type': 'ksponly', + 'ksp_type': 'preonly', + 'pc_type': 'lu', + 'ksp_converged_reason': None, +} +asqstepper = SerialMiniApp(dt, theta, f0, form_mass, form_function, + solver_parameters=params) + +# Run gusto stepper +timestepper.run(0., tmax) +gusto_f_end = timestepper.fields('f') +Print(f"{errornorm(f_end_expr, gusto_f_end) = }") + +# Run asQ stepper +asqstepper.solve(nt=timestepper.step-1) +asq_f_end = asqstepper.w1 +Print(f"{errornorm(f_end_expr, asq_f_end) = }") + +# error +Print(f"{errornorm(gusto_f_end, asq_f_end) = }") diff --git a/case_studies/gusto/gusto_interface.py b/case_studies/gusto/gusto_interface.py new file mode 100644 index 00000000..172f9f51 --- /dev/null +++ b/case_studies/gusto/gusto_interface.py @@ -0,0 +1,66 @@ +from firedrake.fml import ( + Term, all_terms, drop, + replace_subject, replace_test_function +) +from firedrake.formmanipulation import split_form +from gusto.labels import time_derivative, prognostic + + +def asq_forms(equation): + NullTerm = Term(None) + + V = equation.function_space + ncpts = len(V) + residual = equation.residual + + # split equation into mass matrix and linear operator + mass = residual.label_map( + lambda t: t.has_label(time_derivative), + map_if_false=drop) + + function = residual.label_map( + lambda t: t.has_label(time_derivative), + map_if_true=drop) + + # generate ufl for mass matrix over given trial/tests + def form_mass(*trials_and_tests): + trials = trials_and_tests[:ncpts] if ncpts > 1 else trials_and_tests[0] + tests = trials_and_tests[ncpts:] if ncpts > 1 else trials_and_tests[1] + + m = mass.label_map( + all_terms, + replace_test_function(tests)) + m = m.label_map( + all_terms, + replace_subject(trials)) + return m.form + + # generate ufl for linear operator over given trial/tests + def form_function(*trials_and_tests): + trials = trials_and_tests[:ncpts] if ncpts > 1 else trials_and_tests[0] + tests = trials_and_tests[ncpts:] + + fields = equation.field_names if ncpts > 1 else [equation.field_name] + + f = NullTerm + for i in range(ncpts): + fi = function.label_map( + lambda t: t.get(prognostic) == fields[i], + lambda t: Term( + split_form(t.form)[i].form, + t.labels), + map_if_false=drop) + + fi = fi.label_map( + all_terms, + replace_test_function(tests[i])) + + fi = fi.label_map( + all_terms, + replace_subject(trials)) + + f += fi + f = f.label_map(lambda t: t is NullTerm, drop) + return f.form + + return form_mass, form_function diff --git a/case_studies/gusto/lswe.py b/case_studies/gusto/lswe.py new file mode 100644 index 00000000..0789cefe --- /dev/null +++ b/case_studies/gusto/lswe.py @@ -0,0 +1,92 @@ +from firedrake import ( + PeriodicSquareMesh, SpatialCoordinate, Constant, + sin, cos, pi, as_vector, errornorm, norm) +from gusto import * +from firedrake.petsc import PETSc +Print = PETSc.Sys.Print + +from utils.serial import SerialMiniApp +from gusto_interface import asq_forms + + +dt = 0.001 +tmax = 30*dt +H = 1 +wx = 2 +wy = 1 +g = 1 + +# Domain +mesh_name = 'linear_sw_mesh' +L = 1 +nx = ny = 20 +mesh = PeriodicSquareMesh(nx, ny, L, direction='both', name=mesh_name) +x, y = SpatialCoordinate(mesh) +domain = Domain(mesh, dt, 'BDM', 1) + +# Equation +parameters = ShallowWaterParameters(H=H, g=g) +fexpr = Constant(1) +eqns = LinearShallowWaterEquations(domain, parameters, fexpr=fexpr) + +# I/O +output = OutputParameters( + dirname=str('output')+"/linear_sw_wave", + dumpfreq=1, + checkpoint=True +) +io = IO(domain, output) +transport_methods = [DefaultTransport(eqns, "D")] + +# Timestepper +stepper = Timestepper(eqns, TrapeziumRule(domain), io, transport_methods) + +# ---------------------------------------------------------------------- # +# Initial conditions +# ---------------------------------------------------------------------- # + +eta = sin(2*pi*(x-L/2)*wx)*cos(2*pi*(y-L/2)*wy) - (1/5)*cos(2*pi*(x-L/2)*wx)*sin(4*pi*(y-L/2)*wy) +Dexpr = H + eta + +u = cos(4*pi*(x-L/2)*wx)*cos(2*pi*(y-L/2)*wy) +v = cos(2*pi*(x-L/2)*wx)*cos(4*pi*(y-L/2)*wy) +uexpr = as_vector([u, v]) + +u0 = stepper.fields("u") +D0 = stepper.fields("D") + +u0.project(uexpr) +D0.interpolate(Dexpr) + +# --------------------------------------------------------------------- # +# Run +# --------------------------------------------------------------------- # + +stepper.run(t=0, tmax=tmax) + +gusto_u = stepper.fields('u') +gusto_D = stepper.fields('D') + + +# --------------------------------------------------------------------- # +# Setup asQ stepper +# --------------------------------------------------------------------- # +form_mass, form_function = asq_forms(eqns) +theta = 0.5 +w0 = Function(eqns.function_space) +w0.subfunctions[0].project(uexpr) +w0.subfunctions[1].interpolate(Dexpr) +params = { + 'snes_type': 'ksponly', + 'ksp_type': 'preonly', + 'pc_type': 'lu', + 'ksp_converged_reason': None, +} +asqstepper = SerialMiniApp(dt, theta, w0, form_mass, form_function, + solver_parameters=params) +asqstepper.solve(nt=stepper.step-1) +asq_u = asqstepper.w1.subfunctions[0] +asq_D = asqstepper.w1.subfunctions[1] + +Print(f"{errornorm(asq_u, gusto_u)/norm(gusto_u) = }") +Print(f"{errornorm(asq_D, gusto_D)/norm(gusto_D) = }") From a075e5de315d645c6be0d9ffb4babf43e7b87d0b Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Wed, 24 Jul 2024 16:59:10 +0100 Subject: [PATCH 2/4] Replace the transporting velocity correctly from Gusto equations and test with Williamson5 example --- case_studies/gusto/gusto_interface.py | 14 ++- case_studies/gusto/lswe.py | 4 + case_studies/gusto/swe.py | 158 ++++++++++++++++++++++++++ 3 files changed, 172 insertions(+), 4 deletions(-) create mode 100644 case_studies/gusto/swe.py diff --git a/case_studies/gusto/gusto_interface.py b/case_studies/gusto/gusto_interface.py index 172f9f51..69079829 100644 --- a/case_studies/gusto/gusto_interface.py +++ b/case_studies/gusto/gusto_interface.py @@ -3,12 +3,12 @@ replace_subject, replace_test_function ) from firedrake.formmanipulation import split_form -from gusto.labels import time_derivative, prognostic +from gusto import time_derivative, prognostic, transporting_velocity +from ufl import replace -def asq_forms(equation): +def asq_forms(equation, transport_velocity_index=None): NullTerm = Term(None) - V = equation.function_space ncpts = len(V) residual = equation.residual @@ -38,7 +38,7 @@ def form_mass(*trials_and_tests): # generate ufl for linear operator over given trial/tests def form_function(*trials_and_tests): trials = trials_and_tests[:ncpts] if ncpts > 1 else trials_and_tests[0] - tests = trials_and_tests[ncpts:] + tests = trials_and_tests[ncpts:2*ncpts] fields = equation.field_names if ncpts > 1 else [equation.field_name] @@ -59,6 +59,12 @@ def form_function(*trials_and_tests): all_terms, replace_subject(trials)) + if transport_velocity_index is not None: + transport_velocity = trials[transport_velocity_index] + fi = fi.label_map( + lambda t: t.has_label(transporting_velocity), + map_if_true=lambda t: Term(replace(t.form, {t.get(transporting_velocity): transport_velocity}), t.labels)) + f += fi f = f.label_map(lambda t: t is NullTerm, drop) return f.form diff --git a/case_studies/gusto/lswe.py b/case_studies/gusto/lswe.py index 0789cefe..076bf2f2 100644 --- a/case_studies/gusto/lswe.py +++ b/case_studies/gusto/lswe.py @@ -88,5 +88,9 @@ asq_u = asqstepper.w1.subfunctions[0] asq_D = asqstepper.w1.subfunctions[1] +# ------------------------------------------------------------------------ # +# Compare results +# ------------------------------------------------------------------------ # + Print(f"{errornorm(asq_u, gusto_u)/norm(gusto_u) = }") Print(f"{errornorm(asq_D, gusto_D)/norm(gusto_D) = }") diff --git a/case_studies/gusto/swe.py b/case_studies/gusto/swe.py new file mode 100644 index 00000000..c83f2efc --- /dev/null +++ b/case_studies/gusto/swe.py @@ -0,0 +1,158 @@ +""" +The Williamson 5 shallow-water test case (flow over topography), solved with a +discretisation of the non-linear shallow-water equations. + +This uses an icosahedral mesh of the sphere, and runs a series of resolutions. +""" + +from gusto import * +from firedrake import (IcosahedralSphereMesh, SpatialCoordinate, + errornorm, norm, as_vector, pi, sqrt, min_value) +from firedrake.petsc import PETSc +from firedrake.output import VTKFile +Print = PETSc.Sys.Print + +from utils.serial import SerialMiniApp +from gusto_interface import asq_forms + +# ---------------------------------------------------------------------------- # +# Test case parameters +# ---------------------------------------------------------------------------- # + +hour = 60.*60. +day = 24.*60.*60. +ref_level = 3 +dt = 1800 +# tmax = 50*day +tmax = day/3 +ndumps = 1 + +# setup shallow water parameters +R = 6371220. +H = 5960. + +# setup input that doesn't change with ref level or dt +parameters = ShallowWaterParameters(H=H) + +# ------------------------------------------------------------------------ # +# Set up model objects +# ------------------------------------------------------------------------ # + +# Domain +mesh = IcosahedralSphereMesh(radius=R, + refinement_level=ref_level, degree=2) +x = SpatialCoordinate(mesh) +domain = Domain(mesh, dt, 'BDM', 1) + +# Equation +Omega = parameters.Omega +fexpr = 2*Omega*x[2]/R +lamda, theta, _ = lonlatr_from_xyz(x[0], x[1], x[2]) +R0 = pi/9. +R0sq = R0**2 +lamda_c = -pi/2. +lsq = (lamda - lamda_c)**2 +theta_c = pi/6. +thsq = (theta - theta_c)**2 +rsq = min_value(R0sq, lsq+thsq) +r = sqrt(rsq) +bexpr = 2000 * (1 - r/R0) +eqns = ShallowWaterEquations(domain, parameters, fexpr=fexpr, bexpr=bexpr) + +# I/O +dirname = "williamson_5_ref%s_dt%s" % (ref_level, dt) +# dumpfreq = int(tmax / (ndumps*dt)) +dumpfreq = 1 +output = OutputParameters( + dirname=dirname, + dumplist_latlon=['D'], + dumpfreq=dumpfreq, +) +diagnostic_fields = [Sum('D', 'topography')] +io = IO(domain, output, diagnostic_fields=diagnostic_fields) + +# Transport schemes +transport_methods = [DGUpwind(eqns, "u"), DGUpwind(eqns, "D")] + +params = { + 'snes_monitor': None, + 'snes_converged_reason': None, + 'snes_rtol': 1e-12, + 'snes_ksp_ew': None, + 'snes_ksp_rtol0': 1e-2, + 'ksp_type': 'fgmres', + 'pc_type': 'lu', + 'ksp_monitor': None, + 'ksp_converged_rate': None, + 'snes_lag_preconditioner': -2, + 'snes_lag_preconditioner_persists': None, +} + +# Time stepper +stepper = Timestepper(eqns, TrapeziumRule(domain, solver_parameters=params), + io, transport_methods) + +# ------------------------------------------------------------------------ # +# Setup Gusto solver +# ------------------------------------------------------------------------ # + +u0 = stepper.fields('u') +D0 = stepper.fields('D') +u_max = 20. # Maximum amplitude of the zonal wind (m/s) +uexpr = as_vector([-u_max*x[1]/R, u_max*x[0]/R, 0.0]) +g = parameters.g +Rsq = R**2 +Dexpr = H - ((R * Omega * u_max + 0.5*u_max**2)*x[2]**2/Rsq)/g - bexpr + +u0.project(uexpr) +D0.interpolate(Dexpr) + +Dbar = Function(D0.function_space()).assign(H) +stepper.set_reference_profiles([('D', Dbar)]) + +# --------------------------------------------------------------------- # +# Setup asQ stepper +# --------------------------------------------------------------------- # +form_mass, form_function = asq_forms(eqns, transport_velocity_index=0) +theta = 0.5 +w0 = Function(eqns.function_space) +w0.subfunctions[0].project(uexpr) +w0.subfunctions[1].interpolate(Dexpr) +asqstepper = SerialMiniApp(dt, theta, w0, form_mass, form_function, + solver_parameters=params) + +# ------------------------------------------------------------------------ # +# Run Gusto solver +# ------------------------------------------------------------------------ # + +stepper.run(t=0, tmax=tmax) + +gusto_u = stepper.fields('u') +gusto_D = stepper.fields('D') + +# ------------------------------------------------------------------------ # +# Run asQ solver +# ------------------------------------------------------------------------ # + +asq_file = VTKFile(f"results/{dirname}/asq_timeseries.pvd") +asq_file.write(*asqstepper.w1.subfunctions) + + +def postprocess(app, step, t): + asq_file.write(*asqstepper.w1.subfunctions) + + +asqstepper.solve(nt=stepper.step-1, postproc=postprocess) +asq_u = asqstepper.w1.subfunctions[0] +asq_D = asqstepper.w1.subfunctions[1] + +# ------------------------------------------------------------------------ # +# Compare results +# ------------------------------------------------------------------------ # + +uinit, Dinit = w0.subfunctions +Print(f"{errornorm(asq_u, uinit)/norm(uinit) = }") +Print(f"{errornorm(asq_D, Dinit)/norm(Dinit) = }") + +Print(f"{errornorm(asq_u, gusto_u)/norm(uinit) = }") +Print(f"{errornorm(asq_D, gusto_D)/norm(uinit) = }") From ae9e6ad1cd651227d0aac380dc237967c9e5ef64 Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Wed, 24 Jul 2024 17:02:54 +0100 Subject: [PATCH 3/4] fstring flake8 --- case_studies/gusto/advection.py | 10 +++++----- case_studies/gusto/diffusion.py | 8 ++++---- case_studies/gusto/swe.py | 8 ++++---- 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/case_studies/gusto/advection.py b/case_studies/gusto/advection.py index 325f2f39..55919f9c 100644 --- a/case_studies/gusto/advection.py +++ b/case_studies/gusto/advection.py @@ -137,8 +137,8 @@ def form_function(*trials_and_tests): # Set up asQ stepper # ------------------------------------------------------------------------ # -Print(f"{equation.function_space = }") -Print(f"{V = }") +Print(f"{equation.function_space=}") +Print(f"{V=}") params = { 'snes_type': 'ksponly', 'ksp_type': 'preonly', @@ -162,10 +162,10 @@ def form_function(*trials_and_tests): # ------------------------------------------------------------------------ # gusto_error = norm(stepper.fields('f') - f_end) / norm(f_end) -Print(f"{gusto_error = }") +Print(f"{gusto_error=}") asq_error = norm(asqstep.w1 - f_end) / norm(f_end) -Print(f"{asq_error = }") +Print(f"{asq_error=}") gusto_asq_error = norm(stepper.fields('f') - asqstep.w1) / norm(f_end) -Print(f"{gusto_asq_error = }") +Print(f"{gusto_asq_error=}") diff --git a/case_studies/gusto/diffusion.py b/case_studies/gusto/diffusion.py index 5d55d15e..90ee4abe 100644 --- a/case_studies/gusto/diffusion.py +++ b/case_studies/gusto/diffusion.py @@ -26,7 +26,7 @@ f_end_expr = (1/(1+4*tmax))*f_init**(1/(1+4*tmax)) V = domain.spaces("DG") -Print(f"{V = }") +Print(f"{V=}") mu = 5. @@ -55,12 +55,12 @@ # Run gusto stepper timestepper.run(0., tmax) gusto_f_end = timestepper.fields('f') -Print(f"{errornorm(f_end_expr, gusto_f_end) = }") +Print(f"{errornorm(f_end_expr, gusto_f_end)=}") # Run asQ stepper asqstepper.solve(nt=timestepper.step-1) asq_f_end = asqstepper.w1 -Print(f"{errornorm(f_end_expr, asq_f_end) = }") +Print(f"{errornorm(f_end_expr, asq_f_end)=}") # error -Print(f"{errornorm(gusto_f_end, asq_f_end) = }") +Print(f"{errornorm(gusto_f_end, asq_f_end)=}") diff --git a/case_studies/gusto/swe.py b/case_studies/gusto/swe.py index c83f2efc..cbd2fbdf 100644 --- a/case_studies/gusto/swe.py +++ b/case_studies/gusto/swe.py @@ -151,8 +151,8 @@ def postprocess(app, step, t): # ------------------------------------------------------------------------ # uinit, Dinit = w0.subfunctions -Print(f"{errornorm(asq_u, uinit)/norm(uinit) = }") -Print(f"{errornorm(asq_D, Dinit)/norm(Dinit) = }") +Print(f"{errornorm(asq_u, uinit)/norm(uinit)=}") +Print(f"{errornorm(asq_D, Dinit)/norm(Dinit)=}") -Print(f"{errornorm(asq_u, gusto_u)/norm(uinit) = }") -Print(f"{errornorm(asq_D, gusto_D)/norm(uinit) = }") +Print(f"{errornorm(asq_u, gusto_u)/norm(uinit)=}") +Print(f"{errornorm(asq_D, gusto_D)/norm(uinit)=}") From 22e12817c1494abf2219b1e9b163515523a8313b Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Thu, 15 Aug 2024 13:42:18 +0100 Subject: [PATCH 4/4] updates to gusto interface --- case_studies/gusto/gusto_interface.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/case_studies/gusto/gusto_interface.py b/case_studies/gusto/gusto_interface.py index 69079829..41b21044 100644 --- a/case_studies/gusto/gusto_interface.py +++ b/case_studies/gusto/gusto_interface.py @@ -1,24 +1,26 @@ +try: + from gusto import time_derivative, prognostic, transporting_velocity +except ModuleNotFoundError as err: + raise type(err)("Gusto must be installed to use the asQ-Gusto interface.") from err + from firedrake.fml import ( Term, all_terms, drop, replace_subject, replace_test_function ) from firedrake.formmanipulation import split_form -from gusto import time_derivative, prognostic, transporting_velocity from ufl import replace def asq_forms(equation, transport_velocity_index=None): NullTerm = Term(None) - V = equation.function_space - ncpts = len(V) - residual = equation.residual + ncpts = len(equation.function_space) # split equation into mass matrix and linear operator - mass = residual.label_map( + mass = equation.residual.label_map( lambda t: t.has_label(time_derivative), map_if_false=drop) - function = residual.label_map( + function = equation.residual.label_map( lambda t: t.has_label(time_derivative), map_if_true=drop)