Skip to content

Commit

Permalink
fix rayleigh friction by adding brackets to interpolator
Browse files Browse the repository at this point in the history
  • Loading branch information
tommbendall committed Dec 18, 2024
1 parent c5f2f9a commit e54806a
Show file tree
Hide file tree
Showing 2 changed files with 10 additions and 10 deletions.
16 changes: 8 additions & 8 deletions gusto/physics/held_suarez_forcing.py
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ def __init__(self, equation, variable_name, parameters, hs_parameters=None):
thermodynamics.exner_pressure(equation.parameters,
self.rho_averaged, self.theta), self.exner)
self.sigma = Function(Vt)
self.kappa = equation.parameters.kappa
kappa = equation.parameters.kappa

T0surf = hs_parameters.T0surf
T0horiz = hs_parameters.T0horiz
Expand All @@ -66,11 +66,11 @@ def __init__(self, equation, variable_name, parameters, hs_parameters=None):
tau_d = hs_parameters.tau_d
tau_u = hs_parameters.tau_u

theta_condition = (T0surf - T0horiz * sin(lat)**2 - (T0vert * ln(self.exner) * cos(lat)**2)/self.kappa)
theta_condition = (T0surf - T0horiz * sin(lat)**2 - (T0vert * ln(self.exner) * cos(lat)**2)/kappa)
Theta_eq = conditional(T0stra/self.exner >= theta_condition, T0stra/self.exner, theta_condition)

# timescale of temperature forcing
tau_cond = (self.sigma**(1/self.kappa) - sigma_b) / (1 - sigma_b)
tau_cond = (self.sigma**(1/kappa) - sigma_b) / (1 - sigma_b)
newton_freq = 1 / tau_d + (1/tau_u - 1/tau_d) * conditional(0 >= tau_cond, 0, tau_cond) * cos(lat)**4
forcing_expr = newton_freq * (self.theta - Theta_eq)

Expand All @@ -80,7 +80,7 @@ def __init__(self, equation, variable_name, parameters, hs_parameters=None):

# Add relaxation term to residual
test = equation.tests[theta_idx]
dx_reduced = dx(degree=4)
dx_reduced = dx(degree=equation.domain.max_quad_degree)
forcing_form = test * self.source_relaxation * dx_reduced
equation.residual += self.label(subject(prognostic(forcing_form, 'theta'), X), self.evaluate)

Expand Down Expand Up @@ -149,10 +149,10 @@ def __init__(self, equation, hs_parameters=None):

self.sigma = Function(Vt)
sigmab = hs_parameters.sigmab
self.kappa = equation.parameters.kappa
kappa = equation.parameters.kappa
tau_fric = 24 * 60 * 60

tau_cond = (self.sigma**(1/self.kappa) - sigmab) / (1 - sigmab)
tau_cond = (self.sigma**(1/kappa) - sigmab) / (1 - sigmab)
wind_timescale = conditional(ge(0, tau_cond), 0, tau_cond) / tau_fric
forcing_expr = u_hori * wind_timescale

Expand All @@ -161,7 +161,7 @@ def __init__(self, equation, hs_parameters=None):

tests = equation.tests
test = tests[u_idx]
dx_reduced = dx(degree=4)
dx_reduced = dx(degree=equation.domain.max_quad_degree)
source_form = inner(test, self.source_friction) * dx_reduced
equation.residual += self.label(subject(prognostic(source_form, 'u'), X), self.evaluate)

Expand All @@ -177,7 +177,7 @@ def evaluate(self, x_in, dt):
"""
self.X.assign(x_in)
self.rho_recoverer.project()
self.exner_interpolator.interpolate
self.exner_interpolator.interpolate()
# Determine sigma:= exner / exner_surf
exner_columnwise, index_data = self.domain.coords.get_column_data(self.exner, self.domain)
sigma_columnwise = np.zeros_like(exner_columnwise)
Expand Down
4 changes: 2 additions & 2 deletions integration-tests/physics/test_held_suarez_friction.py
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ def run_apply_rayleigh_friction(dirname):

# Answer: slower winds than initially
u_true = Function(Vu)
u_true.project(as_vector([(864 - 864/86400), 0.0]))
u_true.project(as_vector([828, 0.0]))

# ------------------------------------------------------------------------ #
# Run
Expand Down Expand Up @@ -110,5 +110,5 @@ def test_rayleigh_friction(tmpdir):
u_z_true = Function(DG0).project(dot(u_true, e_z))

denom = norm(u_x_true)
assert norm(u_x_final - u_x_true) / denom < 0.01, 'Final horizontal wind is incorrect'
assert norm(u_x_final - u_x_true) / denom < 0.0001, 'Final horizontal wind is incorrect'
assert norm(u_z_final - u_z_true) < 1e-12, 'Final vertical wind is incorrect'

0 comments on commit e54806a

Please sign in to comment.