From e433e7796a9bab34ab5d7bb2002a79a953da182f Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Tue, 21 Nov 2023 22:49:30 +0000 Subject: [PATCH 1/9] SplitPrescribedTransport test of transporting velocity --- .../model/test_prescribed_transport.py | 25 +++++++++++++------ 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/integration-tests/model/test_prescribed_transport.py b/integration-tests/model/test_prescribed_transport.py index cf31b2bfc..f63c0dc61 100644 --- a/integration-tests/model/test_prescribed_transport.py +++ b/integration-tests/model/test_prescribed_transport.py @@ -1,19 +1,21 @@ """ -This tests the prescribed wind feature of the PrescribedTransport timestepper. +This tests the prescribed wind feature of the PrescribedTransport and +SplitPresribedTransport (with no physics schemes) timesteppers. A tracer is transported with a time-varying wind and the computed solution is compared with a true one to check that the transport is working correctly. """ from gusto import * from firedrake import sin, cos, norm, pi, as_vector +import pytest def run(timestepper, tmax, f_end): timestepper.run(0, tmax) return norm(timestepper.fields("f") - f_end) / norm(f_end) - -def test_prescribed_transport_setup(tmpdir, tracer_setup): +@pytest.mark.parametrize('timestep_method', ['prescribed', 'split_prescribed']) +def test_prescribed_transport_setup(tmpdir, tracer_setup, timestep_method): # Make domain using routine from conftest geometry = "slice" @@ -31,11 +33,20 @@ def u_evaluation(t): sin(2*pi*t/setup.tmax)*sin(pi*z)]) transport_scheme = SSPRK3(domain) - transport_method = DGUpwind(eqn, 'f') - timestepper = PrescribedTransport(eqn, transport_scheme, setup.io, - transport_method, - prescribed_transporting_velocity=u_evaluation) + if timestep_method == 'prescribed': + transport_method = DGUpwind(eqn, 'f') + timestepper = PrescribedTransport(eqn, transport_scheme, setup.io, + transport_method, + prescribed_transporting_velocity=u_evaluation) + elif timestep_method == 'split_prescribed': + transport_method = [DGUpwind(eqn, 'f')] + timestepper = SplitPrescribedTransport(eqn, transport_scheme, setup.io, + transport_method, + prescribed_transporting_velocity=u_evaluation) + else: + raise NotImplementedError + # Initial conditions timestepper.fields("f").interpolate(setup.f_init) From 730b1074a1c9f9ac6d1176e766a03cf84c98d4a4 Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Tue, 21 Nov 2023 23:01:04 +0000 Subject: [PATCH 2/9] preciptation integration test now uses SplitPrescribedTransport --- integration-tests/model/test_prescribed_transport.py | 4 ++-- integration-tests/physics/test_precipitation.py | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/integration-tests/model/test_prescribed_transport.py b/integration-tests/model/test_prescribed_transport.py index f63c0dc61..30cbc0292 100644 --- a/integration-tests/model/test_prescribed_transport.py +++ b/integration-tests/model/test_prescribed_transport.py @@ -1,5 +1,5 @@ """ -This tests the prescribed wind feature of the PrescribedTransport and +This tests the prescribed wind feature of the PrescribedTransport and SplitPresribedTransport (with no physics schemes) timesteppers. A tracer is transported with a time-varying wind and the computed solution is compared with a true one to check that the transport is working correctly. @@ -14,6 +14,7 @@ def run(timestepper, tmax, f_end): timestepper.run(0, tmax) return norm(timestepper.fields("f") - f_end) / norm(f_end) + @pytest.mark.parametrize('timestep_method', ['prescribed', 'split_prescribed']) def test_prescribed_transport_setup(tmpdir, tracer_setup, timestep_method): @@ -47,7 +48,6 @@ def u_evaluation(t): else: raise NotImplementedError - # Initial conditions timestepper.fields("f").interpolate(setup.f_init) timestepper.fields("u").project(u_evaluation(Constant(0.0))) diff --git a/integration-tests/physics/test_precipitation.py b/integration-tests/physics/test_precipitation.py index 4c4a6441b..6f1b340e7 100644 --- a/integration-tests/physics/test_precipitation.py +++ b/integration-tests/physics/test_precipitation.py @@ -47,12 +47,12 @@ def setup_fallout(dirname): # Physics schemes rainfall_method = DGUpwind(eqn, 'rain', outflow=True) - physics_parametrisations = [Fallout(eqn, 'rain', domain, rainfall_method)] + physics_parametrisations = [(Fallout(eqn, 'rain', domain, rainfall_method), SSPRK3(domain))] # build time stepper scheme = SSPRK3(domain) - stepper = PrescribedTransport(eqn, scheme, io, transport_method, - physics_parametrisations=physics_parametrisations) + stepper = SplitPrescribedTransport(eqn, scheme, io, transport_method, + physics_schemes=physics_parametrisations) # ------------------------------------------------------------------------ # # Initial conditions From 4a2df00ecddf6a10f45b88a6a6c63b03a4525b43 Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Wed, 22 Nov 2023 02:31:52 +0000 Subject: [PATCH 3/9] added a test for the terminator physics --- .../physics/test_terminator_toy.py | 111 ++++++++++++++++++ integration-tests/physics/test_wind_drag.py | 2 +- 2 files changed, 112 insertions(+), 1 deletion(-) create mode 100644 integration-tests/physics/test_terminator_toy.py diff --git a/integration-tests/physics/test_terminator_toy.py b/integration-tests/physics/test_terminator_toy.py new file mode 100644 index 000000000..39e85d938 --- /dev/null +++ b/integration-tests/physics/test_terminator_toy.py @@ -0,0 +1,111 @@ +""" +This tests the terminator toy physics scheme that models the interaction +of two chemical species through coupled ODEs. +""" + +from gusto import * +from firedrake import IcosahedralSphereMesh, Constant, ge, le, exp, cos, \ + sin, conditional, interpolate, SpatialCoordinate, VectorFunctionSpace, \ + Function, assemble, dx, FunctionSpace, pi, max_value, acos, as_vector, \ + errornorm, norm +import numpy as np + +def run_terminator_toy(dirname): + + # ------------------------------------------------------------------------ # + # Set up model objects + # ------------------------------------------------------------------------ # + + # A much larger timestep than in proper simulations, as this + # tests moving towards a steady state with no flow. + dt = 100000. + + # Make the mesh and domain + R = 6371220. + mesh = IcosahedralSphereMesh(radius=R, + refinement_level=3, degree=2) + + domain = Domain(mesh, dt, 'BDM', 1) + + # get lat lon coordinates + x = SpatialCoordinate(mesh) + lamda, theta, _ = lonlatr_from_xyz(x[0], x[1], x[2]) + + domain = Domain(mesh, dt, 'BDM', 1) + + # Define the interacting species + X = ActiveTracer(name='X', space='DG', + variable_type=TracerVariableType.mixing_ratio, + transport_eqn=TransportEquationType.advective) + + X2 = ActiveTracer(name='X2', space='DG', + variable_type=TracerVariableType.mixing_ratio, + transport_eqn=TransportEquationType.advective) + + tracers = [X, X2] + + # Equation + V = domain.spaces("HDiv") + eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu = V) + + output = OutputParameters(dirname=dirname+"/terminator_toy", + dumpfreq=10) + io = IO(domain, output) + + # Define the reaction rates: + theta_c = np.pi/9. + lamda_c = -np.pi/3. + + k1 = max_value(0, sin(theta)*sin(theta_c) + cos(theta)*cos(theta_c)*cos(lamda-lamda_c)) + k2 = 1 + + physics_schemes = [(TerminatorToy(eqn, k1=k1, k2=k2, species1_name='X', + species2_name='X2'), BackwardEuler(domain))] + + # Set up a non-divergent, time-varying, velocity field + def u_t(t): + return as_vector((Constant(0)*lamda, Constant(0)*lamda, Constant(0)*lamda)) + + transport_scheme = SSPRK3(domain) + transport_method = [DGUpwind(eqn, 'X'), DGUpwind(eqn, 'X2')] + + X_0 = 4e-6 + 0*lamda + X2_0 = 0*lamda + + transport_scheme = SSPRK3(domain)#, limiter=MixedLimiter) + transport_method = [DGUpwind(eqn, 'X'), DGUpwind(eqn, 'X2')] + + stepper = SplitPrescribedTransport(eqn, transport_scheme, io, + spatial_methods=transport_method, + physics_schemes=physics_schemes, + prescribed_transporting_velocity=u_t) + + stepper.fields("X").interpolate(X_0) + stepper.fields("X2").interpolate(X2_0) + + #True fields to compare to + + stepper.run(t=0, tmax=10*dt) + + steady_space = domain.spaces('DG') + X_steady = Function(steady_space) + X2_steady = Function(steady_space) + + X_T_0 = 4e-6 + r = k1/(4*k2) + D_val = sqrt(r**2 + 2*X_T_0*r) + + X_steady.interpolate(D_val - r) + X2_steady.interpolate(0.5*(X_T_0 - D_val + r)) + + return stepper, X_steady, X2_steady + +def test_terminator_toy_setup(tmpdir): + dirname=str(tmpdir) + stepper, X_steady, X2_steady = run_terminator_toy(dirname) + X_field = stepper.fields("X") + X2_field = stepper.fields("X2") + + #Compute + assert errornorm(X_field, X_steady)/norm(X_steady) < 0.25, "The X field is not sufficiently close to the steady state profile" + assert errornorm(X2_field, X2_steady)/norm(X2_steady) < 0.25, "The X2 field is not sufficiently close to the steady state profile" diff --git a/integration-tests/physics/test_wind_drag.py b/integration-tests/physics/test_wind_drag.py index 6bf107a9d..c6aec9bc4 100644 --- a/integration-tests/physics/test_wind_drag.py +++ b/integration-tests/physics/test_wind_drag.py @@ -37,7 +37,7 @@ def run_wind_drag(dirname, implicit_formulation): eqn = CompressibleEulerEquations(domain, parameters) # I/O - output = OutputParameters(dirname=dirname+"/surface_fluxes", + output = OutputParameters(dirname=dirname+"/wind_drag", dumpfreq=1, dumplist=['u']) io = IO(domain, output) From 3076528e8720b4b4a8265260d5d614dc84263d05 Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Wed, 22 Nov 2023 22:30:51 +0000 Subject: [PATCH 4/9] lint terminator physics test --- .../physics/test_terminator_toy.py | 71 ++++++++++--------- 1 file changed, 36 insertions(+), 35 deletions(-) diff --git a/integration-tests/physics/test_terminator_toy.py b/integration-tests/physics/test_terminator_toy.py index 39e85d938..5521bdef2 100644 --- a/integration-tests/physics/test_terminator_toy.py +++ b/integration-tests/physics/test_terminator_toy.py @@ -4,19 +4,19 @@ """ from gusto import * -from firedrake import IcosahedralSphereMesh, Constant, ge, le, exp, cos, \ - sin, conditional, interpolate, SpatialCoordinate, VectorFunctionSpace, \ - Function, assemble, dx, FunctionSpace, pi, max_value, acos, as_vector, \ +from firedrake import IcosahedralSphereMesh, Constant, cos, \ + sin, SpatialCoordinate, Function, max_value, as_vector, \ errornorm, norm import numpy as np + def run_terminator_toy(dirname): # ------------------------------------------------------------------------ # # Set up model objects # ------------------------------------------------------------------------ # - # A much larger timestep than in proper simulations, as this + # A much larger timestep than in proper simulations, as this # tests moving towards a steady state with no flow. dt = 100000. @@ -24,88 +24,89 @@ def run_terminator_toy(dirname): R = 6371220. mesh = IcosahedralSphereMesh(radius=R, refinement_level=3, degree=2) - + domain = Domain(mesh, dt, 'BDM', 1) - + # get lat lon coordinates x = SpatialCoordinate(mesh) lamda, theta, _ = lonlatr_from_xyz(x[0], x[1], x[2]) domain = Domain(mesh, dt, 'BDM', 1) - + # Define the interacting species X = ActiveTracer(name='X', space='DG', - variable_type=TracerVariableType.mixing_ratio, - transport_eqn=TransportEquationType.advective) + variable_type=TracerVariableType.mixing_ratio, + transport_eqn=TransportEquationType.advective) X2 = ActiveTracer(name='X2', space='DG', - variable_type=TracerVariableType.mixing_ratio, - transport_eqn=TransportEquationType.advective) - + variable_type=TracerVariableType.mixing_ratio, + transport_eqn=TransportEquationType.advective) + tracers = [X, X2] - + # Equation V = domain.spaces("HDiv") - eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu = V) - + eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu=V) + output = OutputParameters(dirname=dirname+"/terminator_toy", dumpfreq=10) io = IO(domain, output) - + # Define the reaction rates: theta_c = np.pi/9. lamda_c = -np.pi/3. - + k1 = max_value(0, sin(theta)*sin(theta_c) + cos(theta)*cos(theta_c)*cos(lamda-lamda_c)) k2 = 1 - + physics_schemes = [(TerminatorToy(eqn, k1=k1, k2=k2, species1_name='X', species2_name='X2'), BackwardEuler(domain))] - + # Set up a non-divergent, time-varying, velocity field def u_t(t): return as_vector((Constant(0)*lamda, Constant(0)*lamda, Constant(0)*lamda)) - + transport_scheme = SSPRK3(domain) transport_method = [DGUpwind(eqn, 'X'), DGUpwind(eqn, 'X2')] - + X_0 = 4e-6 + 0*lamda X2_0 = 0*lamda - - transport_scheme = SSPRK3(domain)#, limiter=MixedLimiter) + + transport_scheme = SSPRK3(domain) transport_method = [DGUpwind(eqn, 'X'), DGUpwind(eqn, 'X2')] - stepper = SplitPrescribedTransport(eqn, transport_scheme, io, + stepper = SplitPrescribedTransport(eqn, transport_scheme, io, spatial_methods=transport_method, physics_schemes=physics_schemes, prescribed_transporting_velocity=u_t) - + stepper.fields("X").interpolate(X_0) stepper.fields("X2").interpolate(X2_0) - - #True fields to compare to - + stepper.run(t=0, tmax=10*dt) - + + # Compute the steady state solution to compare to steady_space = domain.spaces('DG') X_steady = Function(steady_space) X2_steady = Function(steady_space) - + X_T_0 = 4e-6 r = k1/(4*k2) D_val = sqrt(r**2 + 2*X_T_0*r) - + X_steady.interpolate(D_val - r) X2_steady.interpolate(0.5*(X_T_0 - D_val + r)) - + return stepper, X_steady, X2_steady - + + def test_terminator_toy_setup(tmpdir): dirname=str(tmpdir) stepper, X_steady, X2_steady = run_terminator_toy(dirname) X_field = stepper.fields("X") X2_field = stepper.fields("X2") - - #Compute + + # Assert that the physics scheme has sufficiently moved + # the species fields near their steady state solutions assert errornorm(X_field, X_steady)/norm(X_steady) < 0.25, "The X field is not sufficiently close to the steady state profile" assert errornorm(X2_field, X2_steady)/norm(X2_steady) < 0.25, "The X2 field is not sufficiently close to the steady state profile" From 2918b52b39b34bbca8ddd136e9106f851e7aeb6b Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Wed, 22 Nov 2023 22:32:36 +0000 Subject: [PATCH 5/9] lint --- integration-tests/physics/test_terminator_toy.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/integration-tests/physics/test_terminator_toy.py b/integration-tests/physics/test_terminator_toy.py index 5521bdef2..77660a8d5 100644 --- a/integration-tests/physics/test_terminator_toy.py +++ b/integration-tests/physics/test_terminator_toy.py @@ -101,12 +101,12 @@ def u_t(t): def test_terminator_toy_setup(tmpdir): - dirname=str(tmpdir) + dirname = str(tmpdir) stepper, X_steady, X2_steady = run_terminator_toy(dirname) X_field = stepper.fields("X") X2_field = stepper.fields("X2") - # Assert that the physics scheme has sufficiently moved - # the species fields near their steady state solutions + # Assert that the physics scheme has sufficiently moved + # the species fields near their steady state solutions assert errornorm(X_field, X_steady)/norm(X_steady) < 0.25, "The X field is not sufficiently close to the steady state profile" assert errornorm(X2_field, X2_steady)/norm(X2_steady) < 0.25, "The X2 field is not sufficiently close to the steady state profile" From 78a0c0fae8ecd6a7d3166c1f0fca38f582b02219 Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Thu, 7 Dec 2023 00:09:29 +0000 Subject: [PATCH 6/9] ref_level=2 for terminator test --- .../physics/test_terminator_toy.py | 22 +++++++++---------- 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/integration-tests/physics/test_terminator_toy.py b/integration-tests/physics/test_terminator_toy.py index 77660a8d5..d09cc1644 100644 --- a/integration-tests/physics/test_terminator_toy.py +++ b/integration-tests/physics/test_terminator_toy.py @@ -18,14 +18,12 @@ def run_terminator_toy(dirname): # A much larger timestep than in proper simulations, as this # tests moving towards a steady state with no flow. - dt = 100000. + dt = 50000. # Make the mesh and domain R = 6371220. mesh = IcosahedralSphereMesh(radius=R, - refinement_level=3, degree=2) - - domain = Domain(mesh, dt, 'BDM', 1) + refinement_level=2, degree=2) # get lat lon coordinates x = SpatialCoordinate(mesh) @@ -64,12 +62,10 @@ def run_terminator_toy(dirname): # Set up a non-divergent, time-varying, velocity field def u_t(t): - return as_vector((Constant(0)*lamda, Constant(0)*lamda, Constant(0)*lamda)) - - transport_scheme = SSPRK3(domain) - transport_method = [DGUpwind(eqn, 'X'), DGUpwind(eqn, 'X2')] + return as_vector((Constant(0)*lamda, Constant(0)*lamda, Constant(0)*lamda)) - X_0 = 4e-6 + 0*lamda + X_T_0 = 4e-6 + X_0 = X_T_0 + 0*lamda X2_0 = 0*lamda transport_scheme = SSPRK3(domain) @@ -90,7 +86,6 @@ def u_t(t): X_steady = Function(steady_space) X2_steady = Function(steady_space) - X_T_0 = 4e-6 r = k1/(4*k2) D_val = sqrt(r**2 + 2*X_T_0*r) @@ -106,7 +101,10 @@ def test_terminator_toy_setup(tmpdir): X_field = stepper.fields("X") X2_field = stepper.fields("X2") + print(errornorm(X_field, X_steady)/norm(X_steady)) + print(errornorm(X2_field, X2_steady)/norm(X2_steady)) + # Assert that the physics scheme has sufficiently moved # the species fields near their steady state solutions - assert errornorm(X_field, X_steady)/norm(X_steady) < 0.25, "The X field is not sufficiently close to the steady state profile" - assert errornorm(X2_field, X2_steady)/norm(X2_steady) < 0.25, "The X2 field is not sufficiently close to the steady state profile" + assert errornorm(X_field, X_steady)/norm(X_steady) < 0.4, "The X field is not sufficiently close to the steady state profile" + assert errornorm(X2_field, X2_steady)/norm(X2_steady) < 0.4, "The X2 field is not sufficiently close to the steady state profile" From 216c7d973bb4f39c37647fbc4108980837d597ea Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Thu, 7 Dec 2023 00:26:13 +0000 Subject: [PATCH 7/9] lint --- integration-tests/physics/test_terminator_toy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/integration-tests/physics/test_terminator_toy.py b/integration-tests/physics/test_terminator_toy.py index d09cc1644..87eefaa5b 100644 --- a/integration-tests/physics/test_terminator_toy.py +++ b/integration-tests/physics/test_terminator_toy.py @@ -62,7 +62,7 @@ def run_terminator_toy(dirname): # Set up a non-divergent, time-varying, velocity field def u_t(t): - return as_vector((Constant(0)*lamda, Constant(0)*lamda, Constant(0)*lamda)) + return as_vector([Constant(0)*lamda,Constant(0)*lamda,Constant(0)*lamda]) X_T_0 = 4e-6 X_0 = X_T_0 + 0*lamda From 7275a70beffad2d03a9883b0b1697a5e4a7eb3cb Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Thu, 7 Dec 2023 00:27:27 +0000 Subject: [PATCH 8/9] lint --- .../physics/test_terminator_toy.py | 110 ------------------ 1 file changed, 110 deletions(-) diff --git a/integration-tests/physics/test_terminator_toy.py b/integration-tests/physics/test_terminator_toy.py index 87eefaa5b..e69de29bb 100644 --- a/integration-tests/physics/test_terminator_toy.py +++ b/integration-tests/physics/test_terminator_toy.py @@ -1,110 +0,0 @@ -""" -This tests the terminator toy physics scheme that models the interaction -of two chemical species through coupled ODEs. -""" - -from gusto import * -from firedrake import IcosahedralSphereMesh, Constant, cos, \ - sin, SpatialCoordinate, Function, max_value, as_vector, \ - errornorm, norm -import numpy as np - - -def run_terminator_toy(dirname): - - # ------------------------------------------------------------------------ # - # Set up model objects - # ------------------------------------------------------------------------ # - - # A much larger timestep than in proper simulations, as this - # tests moving towards a steady state with no flow. - dt = 50000. - - # Make the mesh and domain - R = 6371220. - mesh = IcosahedralSphereMesh(radius=R, - refinement_level=2, degree=2) - - # get lat lon coordinates - x = SpatialCoordinate(mesh) - lamda, theta, _ = lonlatr_from_xyz(x[0], x[1], x[2]) - - domain = Domain(mesh, dt, 'BDM', 1) - - # Define the interacting species - X = ActiveTracer(name='X', space='DG', - variable_type=TracerVariableType.mixing_ratio, - transport_eqn=TransportEquationType.advective) - - X2 = ActiveTracer(name='X2', space='DG', - variable_type=TracerVariableType.mixing_ratio, - transport_eqn=TransportEquationType.advective) - - tracers = [X, X2] - - # Equation - V = domain.spaces("HDiv") - eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu=V) - - output = OutputParameters(dirname=dirname+"/terminator_toy", - dumpfreq=10) - io = IO(domain, output) - - # Define the reaction rates: - theta_c = np.pi/9. - lamda_c = -np.pi/3. - - k1 = max_value(0, sin(theta)*sin(theta_c) + cos(theta)*cos(theta_c)*cos(lamda-lamda_c)) - k2 = 1 - - physics_schemes = [(TerminatorToy(eqn, k1=k1, k2=k2, species1_name='X', - species2_name='X2'), BackwardEuler(domain))] - - # Set up a non-divergent, time-varying, velocity field - def u_t(t): - return as_vector([Constant(0)*lamda,Constant(0)*lamda,Constant(0)*lamda]) - - X_T_0 = 4e-6 - X_0 = X_T_0 + 0*lamda - X2_0 = 0*lamda - - transport_scheme = SSPRK3(domain) - transport_method = [DGUpwind(eqn, 'X'), DGUpwind(eqn, 'X2')] - - stepper = SplitPrescribedTransport(eqn, transport_scheme, io, - spatial_methods=transport_method, - physics_schemes=physics_schemes, - prescribed_transporting_velocity=u_t) - - stepper.fields("X").interpolate(X_0) - stepper.fields("X2").interpolate(X2_0) - - stepper.run(t=0, tmax=10*dt) - - # Compute the steady state solution to compare to - steady_space = domain.spaces('DG') - X_steady = Function(steady_space) - X2_steady = Function(steady_space) - - r = k1/(4*k2) - D_val = sqrt(r**2 + 2*X_T_0*r) - - X_steady.interpolate(D_val - r) - X2_steady.interpolate(0.5*(X_T_0 - D_val + r)) - - return stepper, X_steady, X2_steady - - -def test_terminator_toy_setup(tmpdir): - dirname = str(tmpdir) - stepper, X_steady, X2_steady = run_terminator_toy(dirname) - X_field = stepper.fields("X") - X2_field = stepper.fields("X2") - - print(errornorm(X_field, X_steady)/norm(X_steady)) - print(errornorm(X2_field, X2_steady)/norm(X2_steady)) - - # Assert that the physics scheme has sufficiently moved - # the species fields near their steady state solutions - assert errornorm(X_field, X_steady)/norm(X_steady) < 0.4, "The X field is not sufficiently close to the steady state profile" - assert errornorm(X2_field, X2_steady)/norm(X2_steady) < 0.4, "The X2 field is not sufficiently close to the steady state profile" From 6200cc0b3a43432202d38f027b994dd13a032b14 Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Tue, 12 Dec 2023 00:58:19 +0000 Subject: [PATCH 9/9] undo file deletion --- .../physics/test_terminator_toy.py | 110 ++++++++++++++++++ 1 file changed, 110 insertions(+) diff --git a/integration-tests/physics/test_terminator_toy.py b/integration-tests/physics/test_terminator_toy.py index e69de29bb..f49f87b62 100644 --- a/integration-tests/physics/test_terminator_toy.py +++ b/integration-tests/physics/test_terminator_toy.py @@ -0,0 +1,110 @@ +""" +This tests the terminator toy physics scheme that models the interaction +of two chemical species through coupled ODEs. +""" + +from gusto import * +from firedrake import IcosahedralSphereMesh, Constant, cos, \ + sin, SpatialCoordinate, Function, max_value, as_vector, \ + errornorm, norm +import numpy as np + + +def run_terminator_toy(dirname): + + # ------------------------------------------------------------------------ # + # Set up model objects + # ------------------------------------------------------------------------ # + + # A much larger timestep than in proper simulations, as this + # tests moving towards a steady state with no flow. + dt = 50000. + + # Make the mesh and domain + R = 6371220. + mesh = IcosahedralSphereMesh(radius=R, + refinement_level=2, degree=2) + + # get lat lon coordinates + x = SpatialCoordinate(mesh) + lamda, theta, _ = lonlatr_from_xyz(x[0], x[1], x[2]) + + domain = Domain(mesh, dt, 'BDM', 1) + + # Define the interacting species + X = ActiveTracer(name='X', space='DG', + variable_type=TracerVariableType.mixing_ratio, + transport_eqn=TransportEquationType.advective) + + X2 = ActiveTracer(name='X2', space='DG', + variable_type=TracerVariableType.mixing_ratio, + transport_eqn=TransportEquationType.advective) + + tracers = [X, X2] + + # Equation + V = domain.spaces("HDiv") + eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu=V) + + output = OutputParameters(dirname=dirname+"/terminator_toy", + dumpfreq=10) + io = IO(domain, output) + + # Define the reaction rates: + theta_c = np.pi/9. + lamda_c = -np.pi/3. + + k1 = max_value(0, sin(theta)*sin(theta_c) + cos(theta)*cos(theta_c)*cos(lamda-lamda_c)) + k2 = 1 + + physics_schemes = [(TerminatorToy(eqn, k1=k1, k2=k2, species1_name='X', + species2_name='X2'), BackwardEuler(domain))] + + # Set up a non-divergent, time-varying, velocity field + def u_t(t): + return as_vector([Constant(0)*lamda, Constant(0)*lamda, Constant(0)*lamda]) + + X_T_0 = 4e-6 + X_0 = X_T_0 + 0*lamda + X2_0 = 0*lamda + + transport_scheme = SSPRK3(domain) + transport_method = [DGUpwind(eqn, 'X'), DGUpwind(eqn, 'X2')] + + stepper = SplitPrescribedTransport(eqn, transport_scheme, io, + spatial_methods=transport_method, + physics_schemes=physics_schemes, + prescribed_transporting_velocity=u_t) + + stepper.fields("X").interpolate(X_0) + stepper.fields("X2").interpolate(X2_0) + + stepper.run(t=0, tmax=10*dt) + + # Compute the steady state solution to compare to + steady_space = domain.spaces('DG') + X_steady = Function(steady_space) + X2_steady = Function(steady_space) + + r = k1/(4*k2) + D_val = sqrt(r**2 + 2*X_T_0*r) + + X_steady.interpolate(D_val - r) + X2_steady.interpolate(0.5*(X_T_0 - D_val + r)) + + return stepper, X_steady, X2_steady + + +def test_terminator_toy_setup(tmpdir): + dirname = str(tmpdir) + stepper, X_steady, X2_steady = run_terminator_toy(dirname) + X_field = stepper.fields("X") + X2_field = stepper.fields("X2") + + print(errornorm(X_field, X_steady)/norm(X_steady)) + print(errornorm(X2_field, X2_steady)/norm(X2_steady)) + + # Assert that the physics scheme has sufficiently moved + # the species fields near their steady state solutions + assert errornorm(X_field, X_steady)/norm(X_steady) < 0.4, "The X field is not sufficiently close to the steady state profile" + assert errornorm(X2_field, X2_steady)/norm(X2_steady) < 0.4, "The X2 field is not sufficiently close to the steady state profile"