From 7e90f9b2254c33e2f2d25d9e9126cf65051aa260 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Fri, 16 Aug 2024 13:55:35 +0100 Subject: [PATCH 1/5] initial changes --- gusto/solvers/linear_solvers.py | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/gusto/solvers/linear_solvers.py b/gusto/solvers/linear_solvers.py index f3a610da7..0ee042ad9 100644 --- a/gusto/solvers/linear_solvers.py +++ b/gusto/solvers/linear_solvers.py @@ -294,6 +294,14 @@ def L_tr(f): else: u_mass = inner(w, (u - u_in))*dx + if equations.X.coriolis: + # this becomes much easier once the PR for vorticity goes in due + # to the new way to define coriolis force + omega = Constant(7.292e-5) + Omega = as_vector([0, 0 , omega]) + coriolis_factor = inner(w, cross(2*Omega, u))*dx + else: + coriolis_factor = None eqn = ( # momentum equation u_mass @@ -317,12 +325,19 @@ def L_tr(f): # condition + dl('+')*jump(u, n=n)*(dS_vp + dS_hp) + dl*dot(u, n)*(ds_tbp + ds_vp) - ) + ) # TODO: can we get this term using FML? # contribution of the sponge term if hasattr(self.equations, "mu"): eqn += dt*self.equations.mu*inner(w, k)*inner(u, k)*dx + + if equations.X.coriolis: + # this becomes much easier once the PR for vorticity goes in due + # to the new way to define coriolis force + omega = Constant(7.292e-5) + Omega = as_vector([0, 0 , omega]) + eqn += inner(w, cross(2*Omega, u))*dx aeqn = lhs(eqn) Leqn = rhs(eqn) @@ -654,7 +669,10 @@ def _setup_solver(self): + inner(phi, (D - D_in)) * dx + beta_d * phi * Dbar * div(u) * dx ) - + + if equation.f: # think this is the same as asking if f exists + f = equation.f + eqn += beta_u_ * f * inner(w, domain.perp(u)) * dx aeqn = lhs(eqn) Leqn = rhs(eqn) From fb773107a53ae04c097951b1683f4572541cf98c Mon Sep 17 00:00:00 2001 From: Witt-D Date: Fri, 16 Aug 2024 15:14:59 +0100 Subject: [PATCH 2/5] first draft of solver change with flake8 --- gusto/solvers/linear_solvers.py | 45 +++++++++++++++------------------ 1 file changed, 21 insertions(+), 24 deletions(-) diff --git a/gusto/solvers/linear_solvers.py b/gusto/solvers/linear_solvers.py index 0ee042ad9..1d4b9bf72 100644 --- a/gusto/solvers/linear_solvers.py +++ b/gusto/solvers/linear_solvers.py @@ -9,8 +9,8 @@ split, LinearVariationalProblem, Constant, LinearVariationalSolver, TestFunctions, TrialFunctions, TestFunction, TrialFunction, lhs, rhs, FacetNormal, div, dx, jump, avg, dS, dS_v, dS_h, ds_v, ds_t, ds_b, - ds_tb, inner, action, dot, grad, Function, VectorSpaceBasis, - BrokenElement, FunctionSpace, MixedFunctionSpace, DirichletBC + ds_tb, inner, action, dot, grad, Function, VectorSpaceBasis, cross, + BrokenElement, FunctionSpace, MixedFunctionSpace, DirichletBC, as_vector ) from firedrake.fml import Term, drop from firedrake.petsc import flatten_parameters @@ -294,14 +294,6 @@ def L_tr(f): else: u_mass = inner(w, (u - u_in))*dx - if equations.X.coriolis: - # this becomes much easier once the PR for vorticity goes in due - # to the new way to define coriolis force - omega = Constant(7.292e-5) - Omega = as_vector([0, 0 , omega]) - coriolis_factor = inner(w, cross(2*Omega, u))*dx - else: - coriolis_factor = None eqn = ( # momentum equation u_mass @@ -324,19 +316,15 @@ def L_tr(f): # through the interior facets and weakly impose the no-slip # condition + dl('+')*jump(u, n=n)*(dS_vp + dS_hp) - + dl*dot(u, n)*(ds_tbp + ds_vp) - ) + + dl*dot(u, n)*(ds_tbp + ds_vp)) # TODO: can we get this term using FML? # contribution of the sponge term if hasattr(self.equations, "mu"): eqn += dt*self.equations.mu*inner(w, k)*inner(u, k)*dx - - if equations.X.coriolis: - # this becomes much easier once the PR for vorticity goes in due - # to the new way to define coriolis force - omega = Constant(7.292e-5) - Omega = as_vector([0, 0 , omega]) + + if equations.parameters.coriolis is not None: + Omega = as_vector([0, 0, equations.parameters.Omega]) eqn += inner(w, cross(2*Omega, u))*dx aeqn = lhs(eqn) @@ -521,6 +509,11 @@ def V(u): if hasattr(self.equations, "mu"): eqn += dt*self.equations.mu*inner(w, k)*inner(u, k)*dx + + if equation.parameters.Omega is not None: + Omega = as_vector((0, 0, equation.parameter.Omega)) + eqn += inner(w, cross(2*Omega, u))*dx + aeqn = lhs(eqn) Leqn = rhs(eqn) @@ -667,12 +660,12 @@ def _setup_solver(self): - beta_u * 0.5 * bbar * div(w*(D-Dbar)) * dx + beta_u * 0.5 * jump((D-Dbar)*w, n) * avg(bbar) * dS + inner(phi, (D - D_in)) * dx - + beta_d * phi * Dbar * div(u) * dx - ) - - if equation.f: # think this is the same as asking if f exists - f = equation.f - eqn += beta_u_ * f * inner(w, domain.perp(u)) * dx + + beta_d * phi * Dbar * div(u) * dx) + + if self.prescribed_fields('coriolis'): + f = self.prescribed_fields('coriolis') + eqn += beta_u_ * f * inner(w, equation.domain.perp(u)) * dx + aeqn = lhs(eqn) Leqn = rhs(eqn) @@ -885,6 +878,10 @@ def _setup_solver(self): + beta_d * phi * Dbar * div(u) * dx ) + if self.prescribed_fields('coriolis'): + f = self.prescribed_fields('coriolis') + eqn += beta_u_ * f * inner(w, equation.domain.perp(u)) * dx + aeqn = lhs(eqn) Leqn = rhs(eqn) From e2d4e7c37261d9ee7637a7cc2b5d76738a756a66 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Fri, 16 Aug 2024 16:31:02 +0100 Subject: [PATCH 3/5] fixed corioilis check for SW and reverted bracket change --- gusto/solvers/linear_solvers.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/gusto/solvers/linear_solvers.py b/gusto/solvers/linear_solvers.py index 1d4b9bf72..6f20f5c3b 100644 --- a/gusto/solvers/linear_solvers.py +++ b/gusto/solvers/linear_solvers.py @@ -316,14 +316,15 @@ def L_tr(f): # through the interior facets and weakly impose the no-slip # condition + dl('+')*jump(u, n=n)*(dS_vp + dS_hp) - + dl*dot(u, n)*(ds_tbp + ds_vp)) + + dl*dot(u, n)*(ds_tbp + ds_vp) + ) # TODO: can we get this term using FML? # contribution of the sponge term if hasattr(self.equations, "mu"): eqn += dt*self.equations.mu*inner(w, k)*inner(u, k)*dx - if equations.parameters.coriolis is not None: + if equations.parameters.Omega is not None: Omega = as_vector([0, 0, equations.parameters.Omega]) eqn += inner(w, cross(2*Omega, u))*dx @@ -660,9 +661,10 @@ def _setup_solver(self): - beta_u * 0.5 * bbar * div(w*(D-Dbar)) * dx + beta_u * 0.5 * jump((D-Dbar)*w, n) * avg(bbar) * dS + inner(phi, (D - D_in)) * dx - + beta_d * phi * Dbar * div(u) * dx) + + beta_d * phi * Dbar * div(u) * dx + ) - if self.prescribed_fields('coriolis'): + if 'coriolis' in self.prescribed_fields._field_names: f = self.prescribed_fields('coriolis') eqn += beta_u_ * f * inner(w, equation.domain.perp(u)) * dx @@ -878,7 +880,7 @@ def _setup_solver(self): + beta_d * phi * Dbar * div(u) * dx ) - if self.prescribed_fields('coriolis'): + if 'coriolis' in self.prescribed_fields._field_names: f = self.prescribed_fields('coriolis') eqn += beta_u_ * f * inner(w, equation.domain.perp(u)) * dx From 5b070b982a254e213a5e67f28e5f331d7dacd687 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Fri, 16 Aug 2024 16:40:34 +0100 Subject: [PATCH 4/5] flake8 --- gusto/solvers/linear_solvers.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gusto/solvers/linear_solvers.py b/gusto/solvers/linear_solvers.py index 6f20f5c3b..5d22b7b82 100644 --- a/gusto/solvers/linear_solvers.py +++ b/gusto/solvers/linear_solvers.py @@ -317,7 +317,7 @@ def L_tr(f): # condition + dl('+')*jump(u, n=n)*(dS_vp + dS_hp) + dl*dot(u, n)*(ds_tbp + ds_vp) - ) + ) # TODO: can we get this term using FML? # contribution of the sponge term @@ -662,7 +662,7 @@ def _setup_solver(self): + beta_u * 0.5 * jump((D-Dbar)*w, n) * avg(bbar) * dS + inner(phi, (D - D_in)) * dx + beta_d * phi * Dbar * div(u) * dx - ) + ) if 'coriolis' in self.prescribed_fields._field_names: f = self.prescribed_fields('coriolis') From 8e56dad67bc78f3f6d3e273bd6f20fd1ff17d006 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Mon, 19 Aug 2024 09:38:48 +0100 Subject: [PATCH 5/5] updated the check for the coriolis force in the solvers --- gusto/solvers/linear_solvers.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/gusto/solvers/linear_solvers.py b/gusto/solvers/linear_solvers.py index 5d22b7b82..8214d2ba2 100644 --- a/gusto/solvers/linear_solvers.py +++ b/gusto/solvers/linear_solvers.py @@ -318,7 +318,6 @@ def L_tr(f): + dl('+')*jump(u, n=n)*(dS_vp + dS_hp) + dl*dot(u, n)*(ds_tbp + ds_vp) ) - # TODO: can we get this term using FML? # contribution of the sponge term if hasattr(self.equations, "mu"): @@ -664,8 +663,8 @@ def _setup_solver(self): + beta_d * phi * Dbar * div(u) * dx ) - if 'coriolis' in self.prescribed_fields._field_names: - f = self.prescribed_fields('coriolis') + if 'coriolis' in equation.prescribed_fields._field_names: + f = equation.prescribed_fields('coriolis') eqn += beta_u_ * f * inner(w, equation.domain.perp(u)) * dx aeqn = lhs(eqn) @@ -880,8 +879,8 @@ def _setup_solver(self): + beta_d * phi * Dbar * div(u) * dx ) - if 'coriolis' in self.prescribed_fields._field_names: - f = self.prescribed_fields('coriolis') + if 'coriolis' in equation.prescribed_fields._field_names: + f = equation.prescribed_fields('coriolis') eqn += beta_u_ * f * inner(w, equation.domain.perp(u)) * dx aeqn = lhs(eqn)