Skip to content

Scalar equation

scalar_equation

Scalar terms (e.g. for temperature and salinity transport).

All terms are considered as if they were on the right-hand side of the equation, leading to the following UFL expression returned by the residual method:

\[ (dq)/dt = sum "term.residual()" \]

This sign convention ensures compatibility with Thetis's time integrators. In general, however, we like to think about the terms as they are on the left-hand side. Therefore, in the residual methods below, we first sum the terms in the variable F as if they were on the left-hand side, i.e.

\[ (dq)/dt + F(q) = 0, \]

and then return -F.

advection_term(eq, trial)

Scalar advection term (non-conservative): u \dot \div(q).

Source code in g-adopt/gadopt/scalar_equation.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def advection_term(
    eq: Equation, trial: Argument | ufl.indexed.Indexed | Function
) -> Form:
    r"""Scalar advection term (non-conservative): u \dot \div(q)."""
    advective_velocity_scaling = getattr(eq, "advective_velocity_scaling", 1)
    u = advective_velocity_scaling * eq.u

    if hasattr(eq, "su_nubar"):  # SU(PG) à la Donea & Huerta (2003)
        phi = eq.test + eq.su_nubar / (dot(u, u) + 1e-12) * dot(u, grad(eq.test))

        # The advection term is not integrated by parts so there are no boundary terms
        F = phi * dot(u, grad(trial)) * eq.dx
    else:
        F = -trial * div(eq.test * u) * eq.dx
        F += trial * dot(eq.n, u) * eq.test * eq.ds  # Boundary term in the weak form

        if not (is_continuous(eq.trial_space) and normal_is_continuous(eq.u)):
            # s = 0: u.n(-) < 0 => flow goes from '+' to '-' => '+' is upwind
            # s = 1: u.n(-) > 0 => flow goes from '-' to '+' => '-' is upwind
            s = 0.5 * (sign(dot(avg(u), eq.n("-"))) + 1.0)
            q_up = trial("-") * s + trial("+") * (1 - s)
            F += jump(eq.test * u, eq.n) * q_up * eq.dS

    for bc_id, bc in eq.bcs.items():
        if "q" in bc:
            # On incoming boundaries, where dot(u, n) < 0, replace trial with bc["q"].
            F += eq.test * min_value(dot(u, eq.n), 0) * (bc["q"] - trial) * eq.ds(bc_id)

    return -F

diffusion_term(eq, trial)

Scalar diffusion term \(-nabla * (kappa grad q)\).

Using the symmetric interior penalty method, the weak form becomes

\[ {:( -int_Omega nabla * (kappa grad q) phi dx , = , int_Omega kappa (grad phi) * (grad q) dx ), ( , - , int_(cc"I" uu cc"I"_v) "jump"(phi bb n) * "avg"(kappa grad q) dS - int_(cc"I" uu cc"I"_v) "jump"(q bb n) * "avg"(kappa grad phi) dS ), ( , + , int_(cc"I" uu cc"I"_v) sigma "avg"(kappa) "jump"(q bb n) * "jump"(phi bb n) dS ) :} \]

where σ is a penalty parameter.

Epshteyn, Y., & Rivière, B. (2007). Estimation of penalty parameters for symmetric interior penalty Galerkin methods. Journal of Computational and Applied Mathematics, 206(2), 843-872.

Source code in g-adopt/gadopt/scalar_equation.py 
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def diffusion_term(
    eq: Equation, trial: Argument | ufl.indexed.Indexed | Function
) -> Form:
    r"""Scalar diffusion term $-nabla * (kappa grad q)$.

    Using the symmetric interior penalty method, the weak form becomes

    $$
    {:( -int_Omega nabla * (kappa grad q) phi dx , = , int_Omega kappa (grad phi) * (grad q) dx ),
      ( , - , int_(cc"I" uu cc"I"_v) "jump"(phi bb n) * "avg"(kappa grad q) dS
          -   int_(cc"I" uu cc"I"_v) "jump"(q bb n) * "avg"(kappa grad phi) dS ),
      ( , + , int_(cc"I" uu cc"I"_v) sigma "avg"(kappa) "jump"(q bb n) * "jump"(phi bb n) dS )
    :}
    $$

    where σ is a penalty parameter.

    Epshteyn, Y., & Rivière, B. (2007).
    Estimation of penalty parameters for symmetric interior penalty Galerkin methods.
    Journal of Computational and Applied Mathematics, 206(2), 843-872.
    """
    kappa = eq.diffusivity
    dim = eq.mesh.geometric_dimension()
    diff_tensor = kappa if len(kappa.ufl_shape) == 2 else kappa * Identity(dim)

    reference_for_diffusion = getattr(eq, "reference_for_diffusion", 0)
    q = trial + reference_for_diffusion

    F = inner(grad(eq.test), dot(diff_tensor, grad(q))) * eq.dx

    sigma = interior_penalty_factor(eq, shift=-1)
    if not is_continuous(eq.trial_space):
        sigma_int = sigma * avg(FacetArea(eq.mesh) / CellVolume(eq.mesh))
        F += (
            sigma_int
            * inner(jump(eq.test, eq.n), dot(avg(diff_tensor), jump(q, eq.n)))
            * eq.dS
        )
        F -= inner(avg(dot(diff_tensor, grad(eq.test))), jump(q, eq.n)) * eq.dS
        F -= inner(jump(eq.test, eq.n), avg(dot(diff_tensor, grad(q)))) * eq.dS

    for bc_id, bc in eq.bcs.items():
        if "q" in bc and "flux" in bc:
            raise ValueError("Cannot apply both `q` and `flux` on the same boundary.")

        if "q" in bc:
            jump_q = trial - bc["q"]
            sigma_ext = sigma * FacetArea(eq.mesh) / CellVolume(eq.mesh)
            # Terms below are similar to the above terms for the DG dS integrals.
            F += (
                2
                * sigma_ext
                * eq.test
                * inner(eq.n, dot(diff_tensor, eq.n))
                * jump_q
                * eq.ds(bc_id)
            )
            F -= inner(dot(diff_tensor, grad(eq.test)), eq.n) * jump_q * eq.ds(bc_id)
            F -= inner(eq.test * eq.n, dot(diff_tensor, grad(q))) * eq.ds(bc_id)
        elif "flux" in bc:
            # Here we need only the third term because we assume jump_q = 0
            # (q_ext = trial) and flux = kappa dq/dn = dot(n, dot(diff_tensor, grad(q)).
            F -= eq.test * bc["flux"] * eq.ds(bc_id)

    return -F

source_term(eq, trial)

Scalar source term s_T.

Source code in g-adopt/gadopt/scalar_equation.py 
127
128
129
130
131
def source_term(eq: Equation, trial: Argument | ufl.indexed.Indexed | Function) -> Form:
    r"""Scalar source term `s_T`."""
    F = -dot(eq.test, eq.source) * eq.dx

    return -F

sink_term(eq, trial)

Scalar sink term \alpha_T T.

Source code in g-adopt/gadopt/scalar_equation.py 
134
135
136
137
138
139
def sink_term(eq: Equation, trial: Argument | ufl.indexed.Indexed | Function) -> Form:
    r"""Scalar sink term `\alpha_T T`."""
    # Implement sink term implicitly at current time step.
    F = dot(eq.test, eq.sink_coeff * trial) * eq.dx

    return -F

mass_term(eq, trial)

UFL form for the mass term used in the time discretisation.

Parameters:

Name Type Description Default
eq Equation

G-ADOPT Equation.

 required
trial Argument | Indexed | Function

Firedrake trial function.

 required

Returns:

Type Description
Form

The UFL form associated with the mass term of the equation.
Source code in g-adopt/gadopt/scalar_equation.py 
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def mass_term(eq: Equation, trial: Argument | ufl.indexed.Indexed | Function) -> Form:
    """UFL form for the mass term used in the time discretisation.

    Args:
        eq:
          G-ADOPT Equation.
        trial:
          Firedrake trial function.

    Returns:
        The UFL form associated with the mass term of the equation.

    """
    return dot(eq.test, trial) * eq.dx