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 = {'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 = {'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, solver_parameters_extra=None, timestepper_kwargs=None, su_advection=False)

Bases: SolverConfigurationMixin, 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	IrksomeIntegrator	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	ConfigType | str | None	Dictionary of solver parameters or a string specifying a default configuration provided to PETSc	None
solver_parameters_extra	ConfigType | None	Dictionary of PETSc solver options used to update the default G-ADOPT options	None
timestepper_kwargs	dict[str, Any] | None	Dictionary of additional keyword arguments passed to the timestepper constructor. Useful for parameterised schemes (e.g., {'order': 5} for IrksomeRadauIIA)	None
su_advection	bool	Boolean activating the streamline-upwind stabilisation scheme when using continuous finite elements	False

Source code in g-adopt/gadopt/transport_solver.py

def __init__(
    self,
    solution: Function,
    /,
    delta_t: Constant,
    timestepper: IrksomeIntegrator,
    *,
    solution_old: Function | None = None,
    eq_attrs: dict[str, float] = {},
    bcs: dict[int, dict[str, Number]] = {},
    solver_parameters: ConfigType | str | None = None,
    solver_parameters_extra: ConfigType | None = None,
    timestepper_kwargs: dict[str, Any] | None = None,
    su_advection: bool = False,
) -> None:
    self.solution = solution
    self.delta_t = delta_t
    self.timestepper = timestepper
    self.timestepper_kwargs = timestepper_kwargs or {}
    self.solution_old = solution_old or Function(solution)
    self.eq_attrs = eq_attrs
    self.bcs = bcs
    self.su_advection = su_advection

    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)

    self.set_boundary_conditions()
    self.set_su_nubar()
    self.set_equation()
    self.set_solver_options(solver_parameters, solver_parameters_extra)
    self.setup_solver()

set_boundary_conditions()

Sets up boundary conditions.

Source code in g-adopt/gadopt/transport_solver.py

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

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 not self.su_advection:
        return

    if (u := getattr(self, "u", self.eq_attrs.get("u"))) is None:
        raise ValueError(
            "'u' must be included into `eq_attrs` if `su_advection` 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.
    kappa = self.eq_attrs.get("diffusivity", 0.0)
    Pe = absv(dot(u, J)) / 2 / (kappa + 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

@abc.abstractmethod
def set_equation(self):
    """Sets up the term contributions in the equation."""
    raise NotImplementedError

set_solver_options(solver_preset, solver_extras=None)

Sets PETSc solver parameters.

Source code in g-adopt/gadopt/transport_solver.py

def set_solver_options(
    self,
    solver_preset: ConfigType | str | None,
    solver_extras: ConfigType | None = None,
) -> None:
    """Sets PETSc solver parameters."""
    if isinstance(solver_preset, Mapping):
        self.add_to_solver_config(solver_preset)
        self.add_to_solver_config(solver_extras)
        self.register_update_callback(self.setup_solver)
        return

    if solver_preset is not None:
        match solver_preset:
            case "direct":
                self.add_to_solver_config(direct_energy_solver_parameters)
            case "iterative":
                self.add_to_solver_config(iterative_energy_solver_parameters)
            case _:
                raise ValueError("Solver type must be 'direct' or 'iterative'.")
    elif self.mesh.topological_dimension == 2:
        self.add_to_solver_config(direct_energy_solver_parameters)
    else:
        self.add_to_solver_config(iterative_energy_solver_parameters)

    if DEBUG >= log_level:
        self.add_to_solver_config({"ksp_monitor": None})
    elif INFO >= log_level:
        self.add_to_solver_config({"ksp_converged_reason": None})

    self.add_to_solver_config(solver_extras)
    self.register_update_callback(self.setup_solver)

setup_solver()

Sets up the timestepper using specified parameters.

Source code in g-adopt/gadopt/transport_solver.py

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.timestepper_kwargs,
    )

solver_callback()

Optional instructions to execute right after a solve.

Source code in g-adopt/gadopt/transport_solver.py

def solver_callback(self) -> None:
    """Optional instructions to execute right after a solve."""
    pass

solve(t=None)

Advances solver in time.

Source code in g-adopt/gadopt/transport_solver.py

def solve(self, t: float | None = None) -> None:
    """Advances solver in time."""
    self.ts.advance(t=t)

    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	IrksomeIntegrator	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
timestepper_kwargs		Dictionary of additional keyword arguments passed to the timestepper constructor. Useful for parameterized schemes (e.g., {'order': 5} for IrksomeRadauIIA)	required
su_advection		Boolean activating the streamline-upwind stabilisation scheme when using continuous finite elements	required

Source code in g-adopt/gadopt/transport_solver.py

def __init__(
    self,
    terms: str | list[str],
    solution: Function,
    /,
    delta_t: Constant,
    timestepper: IrksomeIntegrator,
    **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	IrksomeIntegrator	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
timestepper_kwargs		Dictionary of additional keyword arguments passed to the timestepper constructor. Useful for parameterized schemes (e.g., {'order': 5} for IrksomeRadauIIA)	required
su_advection		Boolean activating the streamline-upwind stabilisation scheme when using continuous finite elements	required

Source code in g-adopt/gadopt/transport_solver.py

def __init__(
    self,
    solution: Function,
    u: Function,
    approximation: BaseApproximation,
    /,
    delta_t: Constant,
    timestepper: IrksomeIntegrator,
    **kwargs,
) -> None:
    self.u = u
    self.approximation = approximation

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

    self.T_old = self.solution_old

DiffusiveSmoothingSolver(solution, wavelength, K=1, bcs=None, solver_parameters=None, integration_quad_degree=None, **kwargs)

Bases: GenericTransportSolver

A class to perform diffusive smoothing by inheriting from GenericTransportSolver.

This class provides functionality to solve a diffusion equation for smoothing a scalar function, using clean inheritance from GenericTransportSolver.

Parameters:

Name	Type	Description	Default
solution	Function	The solution function to store the smoothed result.	required
wavelength	Number	The wavelength for diffusion.	required
K	Union[Function, Number]	Diffusion tensor. Defaults to 1.	1
bcs	dict[int, dict[str, int | float]] | None	Boundary conditions to impose on the solution.	None
solver_parameters	dict[str, str | float] | None	Solver parameters for the solver. Defaults to None.	None
integration_quad_degree	int | None	Quadrature degree for integrating the diffusion tensor. If None, defaults to 2p+1 where p is the polynomial degree.	None
**kwargs		Additional keyword arguments to pass to the solver.	{}

Source code in g-adopt/gadopt/transport_solver.py

def __init__(
    self,
    solution: Function,
    wavelength: Number,
    K: Function | Number = 1,
    bcs: dict[int, dict[str, int | float]] | None = None,
    solver_parameters: dict[str, str | float] | None = None,
    integration_quad_degree: int | None = None,
    **kwargs
):
    # Extract function space from solution
    function_space = solution.function_space()

    # Calculate diffusive time step
    dt = self._calculate_diffusive_time_step(function_space, wavelength, K, integration_quad_degree)

    # Initialise the parent GenericTransportSolver
    super().__init__(
        "diffusion",
        solution,
        dt,
        BackwardEuler,
        eq_attrs={"diffusivity": ensure_constant(K)},
        bcs=bcs,
        solver_parameters=solver_parameters,
        **kwargs
    )

action(field)

Apply smoothing action to an input field.

Parameters:

Name	Type	Description	Default
field	Function	The input field to be smoothed.	required

Note

The smoothed result is stored in the solution function passed to the constructor.

Source code in g-adopt/gadopt/transport_solver.py

def action(self, field: Function) -> None:
    """Apply smoothing action to an input field.

    Args:
        field (firedrake.Function): The input field to be smoothed.

    Note:
        The smoothed result is stored in the solution function passed to the constructor.
    """
    # Start with the input field
    self.solution.assign(field)

    # Solve the diffusion equation (inherited from GenericTransportSolver)
    self.solve()