Skip to content

Transport solver

transport_solver

This module provides a minimal solver, independent of any approximation, for generic transport equations, which may include advection, diffusion, sink, and source terms, and a fine-tuned solver class for the energy conservation equation. Users instantiate the GenericTransportSolver or EnergySolver classes by providing the appropriate documented parameters and call the solve method to request a solver update.

iterative_energy_solver_parameters: dict[str, Any] = {'mat_type': 'aij', 'snes_type': 'ksponly', 'ksp_type': 'gmres', 'ksp_rtol': 1e-05, 'pc_type': 'sor'} module-attribute

Default iterative solver parameters for solution of energy equation.

Configured to use the GMRES Krylov scheme with Successive Over Relaxation (SOR) preconditioning. Note that default energy solver parameters can be augmented or adjusted by accessing the solver_parameter dictionary.

Examples:

>>> energy_solver.solver_parameters['ksp_converged_reason'] = None
>>> energy_solver.solver_parameters['ksp_rtol'] = 1e-4
Note

G-ADOPT defaults to iterative solvers in 3-D.

direct_energy_solver_parameters: dict[str, Any] = {'mat_type': 'aij', 'snes_type': 'ksponly', 'ksp_type': 'preonly', 'pc_type': 'lu', 'pc_factor_mat_solver_type': 'mumps'} module-attribute

Default direct solver parameters for solution of energy equation.

Configured to use LU factorisation performed via the MUMPS library.

Note

G-ADOPT defaults to direct solvers in 2-D.

GenericTransportBase(solution, /, delta_t, timestepper, *, solution_old=None, eq_attrs={}, bcs={}, solver_parameters=None, su_diffusivity=None)

Bases: ABC

Base class for advancing a generic transport equation in time.

All combinations of advection, diffusion, sink, and source terms are handled.

Note: The solution field is updated in place.

Parameters:

Name Type Description Default
solution Function

Firedrake function for the field of interest

required
delta_t Constant

Simulation time step

required
timestepper RungeKuttaTimeIntegrator

Runge-Kutta time integrator employing an explicit or implicit numerical scheme

required
solution_old Function | None

Firedrake function holding the solution field at the previous time step

None
eq_attrs dict[str, float]

Dictionary of terms arguments and their values

{}
bcs dict[int, dict[str, Number]]

Dictionary specifying boundary conditions (identifier, type, and value)

{}
solver_parameters dict[str, str | Number] | str | None

Dictionary of solver parameters or a string specifying a default configuration provided to PETSc

None
su_diffusivity float | None

Float activating the streamline-upwind stabilisation scheme and specifying the corresponding diffusivity

None
Source code in g-adopt/gadopt/transport_solver.py
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def __init__(
    self,
    solution: Function,
    /,
    delta_t: Constant,
    timestepper: RungeKuttaTimeIntegrator,
    *,
    solution_old: Function | None = None,
    eq_attrs: dict[str, float] = {},
    bcs: dict[int, dict[str, Number]] = {},
    solver_parameters: dict[str, str | Number] | str | None = None,
    su_diffusivity: float | None = None,
) -> None:
    self.solution = solution
    self.delta_t = delta_t
    self.timestepper = timestepper
    self.solution_old = solution_old or Function(solution)
    self.eq_attrs = eq_attrs
    self.bcs = bcs
    self.solver_parameters = solver_parameters
    self.su_diffusivity = su_diffusivity

    self.solution_space = solution.function_space()
    self.mesh = self.solution_space.mesh()
    self.test = TestFunction(self.solution_space)

    self.continuous_solution = is_continuous(self.solution)

    # Solver object is set up later to permit editing default solver options.
    self._solver_ready = False

    self.set_boundary_conditions()
    self.set_su_nubar()
    self.set_equation()
    self.set_solver_options()

set_boundary_conditions()

Sets up boundary conditions.

Source code in g-adopt/gadopt/transport_solver.py
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def set_boundary_conditions(self) -> None:
    """Sets up boundary conditions."""
    self.strong_bcs = []
    self.weak_bcs = {}

    for bc_id, bc in self.bcs.items():
        weak_bc = {}

        for bc_type, value in bc.items():
            if bc_type == self.strong_bcs_tag:
                if self.continuous_solution:
                    strong_bc = DirichletBC(self.solution_space, value, bc_id)
                    self.strong_bcs.append(strong_bc)
                else:
                    weak_bc["q"] = value
            else:
                weak_bc[bc_type] = value

        self.weak_bcs[bc_id] = weak_bc

set_su_nubar()

Sets up the advection streamline-upwind scheme (Donea & Huerta, 2003).

Columns of Jacobian J are the vectors that span the quad/hex and can be seen as unit vectors scaled with the dx/dy/dz in that direction (assuming physical coordinates x, y, z aligned with local coordinates). Thus u^T J is (dx * u , dy * v). Following (2.44c), Pe = u^T J / 2 kappa, and beta(Pe) is the xibar vector in (2.44a). Finally, we get the artificial diffusion nubar from (2.49).

Donea, J., & Huerta, A. (2003). Finite element methods for flow problems. John Wiley & Sons.

Source code in g-adopt/gadopt/transport_solver.py
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
def set_su_nubar(self) -> None:
    """Sets up the advection streamline-upwind scheme (Donea & Huerta, 2003).

    Columns of Jacobian J are the vectors that span the quad/hex and can be seen as
    unit vectors scaled with the dx/dy/dz in that direction (assuming physical
    coordinates x, y, z aligned with local coordinates).
    Thus u^T J is (dx * u , dy * v). Following (2.44c), Pe = u^T J / 2 kappa, and
    beta(Pe) is the xibar vector in (2.44a). Finally, we get the artificial
    diffusion nubar from (2.49).

    Donea, J., & Huerta, A. (2003).
    Finite element methods for flow problems.
    John Wiley & Sons.
    """
    if self.su_diffusivity is None:
        return

    if (u := getattr(self, "u", self.eq_attrs.get("u"))) is None:
        raise ValueError(
            "'u' must be included into `eq_attrs` if `su_diffusivity` is given."
        )

    if not self.continuous_solution:
        raise TypeError("SU advection requires a continuous function space.")

    log("Using SU advection")

    J = Function(TensorFunctionSpace(self.mesh, "DQ", 1), name="Jacobian")
    J.interpolate(Jacobian(self.mesh))
    # Calculate grid Peclet number. Note the use of a lower bound for diffusivity if
    # a pure advection scenario is considered.
    Pe = absv(dot(u, J)) / 2 / (self.su_diffusivity + 1e-12)
    beta_Pe = as_vector([1 / tanh(Pe_i + 1e-6) - 1 / (Pe_i + 1e-6) for Pe_i in Pe])
    nubar = dot(absv(dot(u, J)), beta_Pe) / 2  # Calculate SU artificial diffusion

    self.eq_attrs["su_nubar"] = nubar

set_equation() abstractmethod

Sets up the term contributions in the equation.

Source code in g-adopt/gadopt/transport_solver.py
194
195
196
197
@abc.abstractmethod
def set_equation(self):
    """Sets up the term contributions in the equation."""
    raise NotImplementedError

set_solver_options()

Sets PETSc solver parameters.

Source code in g-adopt/gadopt/transport_solver.py
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
def set_solver_options(self) -> None:
    """Sets PETSc solver parameters."""
    if isinstance(self.solver_parameters, dict):
        return

    if self.solver_parameters is not None:
        match self.solver_parameters:
            case "direct":
                self.solver_parameters = direct_energy_solver_parameters.copy()
            case "iterative":
                self.solver_parameters = iterative_energy_solver_parameters.copy()
            case _:
                raise ValueError(
                    f"Solver type '{self.solver_parameters}' not implemented."
                )
    elif self.mesh.topological_dimension() == 2:
        self.solver_parameters = direct_energy_solver_parameters.copy()
    else:
        self.solver_parameters = iterative_energy_solver_parameters.copy()

    if DEBUG >= log_level:
        self.solver_parameters["ksp_monitor"] = None
    elif INFO >= log_level:
        self.solver_parameters["ksp_converged_reason"] = None

setup_solver()

Sets up the timestepper using specified parameters.

Source code in g-adopt/gadopt/transport_solver.py
224
225
226
227
228
229
230
231
232
233
234
235
def setup_solver(self) -> None:
    """Sets up the timestepper using specified parameters."""
    self.ts = self.timestepper(
        self.equation,
        self.solution,
        self.delta_t,
        solution_old=self.solution_old,
        solver_parameters=self.solver_parameters,
        strong_bcs=self.strong_bcs,
    )

    self._solver_ready = True

solver_callback()

Optional instructions to execute right after a solve.

Source code in g-adopt/gadopt/transport_solver.py
237
238
239
def solver_callback(self) -> None:
    """Optional instructions to execute right after a solve."""
    pass

solve(t=0, update_forcings=None)

Advances solver in time.

Source code in g-adopt/gadopt/transport_solver.py
241
242
243
244
245
246
247
248
def solve(self, t: Number = 0, update_forcings: Callable | None = None) -> None:
    """Advances solver in time."""
    if not self._solver_ready:
        self.setup_solver()

    self.ts.advance(t, update_forcings)

    self.solver_callback()

GenericTransportSolver(terms, solution, /, delta_t, timestepper, **kwargs)

Bases: GenericTransportBase

Advances in time a generic transport equation.

Note: The solution field is updated in place.

Terms and Attributes

This solver handles all combinations of advection, diffusion, sink, and source terms. Depending on the included terms, specific attributes must be provided according to:

Term Required attribute(s) Optional attribute(s)
advection u advective_velocity_scaling, su_nubar
diffusion diffusivity reference_for_diffusion, interior_penalty
source source
sink sink_coeff

Parameters:

Name Type Description Default
terms str | list[str]

List of equation terms to include (a string for a single term is accepted)

required
solution Function

Firedrake function for the field of interest

required
delta_t Constant

Simulation time step

required
timestepper RungeKuttaTimeIntegrator

Runge-Kutta time integrator employing an explicit or implicit numerical scheme

required
solution_old

Firedrake function holding the solution field at the previous time step

required
eq_attrs

Dictionary of terms arguments and their values

required
bcs

Dictionary specifying boundary conditions (identifier, type, and value)

required
solver_parameters

Dictionary of solver parameters or a string specifying a default configuration provided to PETSc

required
su_diffusivity

Float activating the streamline-upwind stabilisation scheme and specifying the corresponding diffusivity

required
Source code in g-adopt/gadopt/transport_solver.py
294
295
296
297
298
299
300
301
302
303
304
305
def __init__(
    self,
    terms: str | list[str],
    solution: Function,
    /,
    delta_t: Constant,
    timestepper: RungeKuttaTimeIntegrator,
    **kwargs,
) -> None:
    self.terms = [terms] if isinstance(terms, str) else terms

    super().__init__(solution, delta_t, timestepper, **kwargs)

EnergySolver(solution, u, approximation, /, delta_t, timestepper, **kwargs)

Bases: GenericTransportBase

Advances in time the energy conservation equation.

Note: The solution field is updated in place.

Parameters:

Name Type Description Default
solution Function

Firedrake function for temperature

required
u Function

Firedrake function for velocity

required
approximation BaseApproximation

G-ADOPT approximation defining terms in the system of equations

required
delta_t Constant

Simulation time step

required
timestepper RungeKuttaTimeIntegrator

Runge-Kutta time integrator employing an explicit or implicit numerical scheme

required
solution_old

Firedrake function holding the solution field at the previous time step

required
bcs

Dictionary specifying boundary conditions (identifier, type, and value)

required
solver_parameters

Dictionary of solver parameters or a string specifying a default configuration provided to PETSc

required
su_diffusivity

Float activating the streamline-upwind stabilisation scheme and specifying the corresponding diffusivity

required
Source code in g-adopt/gadopt/transport_solver.py
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
def __init__(
    self,
    solution: Function,
    u: Function,
    approximation: BaseApproximation,
    /,
    delta_t: Constant,
    timestepper: RungeKuttaTimeIntegrator,
    **kwargs,
) -> None:
    self.u = u
    self.approximation = approximation

    super().__init__(solution, delta_t, timestepper, **kwargs)

    self.T_old = self.solution_old