Scalar-Transport

Steady Scalar-Transport

class flatiron_tk.physics.steady_scalar_transport.SteadyScalarTransport(mesh, tag='u', q_degree=4, **kwargs)

Bases: PhysicsProblem

A class to represent the steady scalar transport equations. This class extends the PhysicsProblem class to include advection, diffusion, and reaction terms.

add_edge_stab(gamma)

Adds edge stabilization to the weak form.

Parameters:

gamma (The edge stabilization parameter.)

add_stab(tau_type='shakib')

Adds the stabilization term to the weak form.

Parameters:

tau_type (The type of stabilization constant to use. Options are 'shakib','su', or 'codina'.)

build_function_space(*args, **kwargs)

Build the function space for the steady scalar transport problem.

Parameters:

• *args (Arguments to be passed to the parent class method.)

• **kwargs (Keyword arguments to be passed to the parent class method.)

flux(q)

Computes the flux of the problem.

Parameters:

q (The flux vector)

Advection_Matrix + Diffusion_Matrix - Reaction_Matrix - Flux_Term = 0 => we subtract the form

get_advection_velocity()

Get the advection velocity for the steady scalar transport problem.

Returns:

The advection velocity function.

get_diffusivity()

Get the diffusivity for the steady scalar transport problem.

Returns:

The diffusivity function.

get_edge_stab(gamma, c)

Computes the edge stabilization term.

Parameters:

• gamma (The edge stabilization parameter.)

• c (The solution function.)

Returns:

J: The edge stabilization term.

get_reaction()

Get the reaction for the steady scalar transport problem.

Returns:

The reaction function.

get_residual()

Computes the residual of the problem.

Returns:

residual: The residual of the problem as a fenicsx form

set_advection_velocity(u)

Set the advection velocity for the steady scalar transport problem.

Parameters:

u (The advection velocity function.)

set_diffusivity(D)

Set the diffusivity for the steady scalar transport problem.

Parameters:

D (The diffusivity function.)

set_reaction(R)

Set the reaction for the steady scalar transport problem.

Parameters:

R (The reaction function.)

set_weak_form(stab=False)

Set the weak form for the steady scalar transport problem.

Parameters:

stab (Whether to include stabilization terms.)

stabilization_constant(tau_type)

Computes the stabilization constant.

Parameters:

tau_type (The type of stabilization constant to use. Options are 'shakib','su', or 'codina'.)

Returns:

tau: The stabilization constant.

This class defines the scalar transport problem, the so-called Advection-Diffusion-Reaction problem

Strong formulation

Define a total computational domian as \(\Omega\) and the Dirichlet and Neumann boundaries \(\Gamma_D\) and \(\Gamma_N\) where \(\Gamma_D \cup \Gamma_D = \partial\Omega\) and \(\Gamma_D \cap \Gamma_D = 0\). Let the trial function space be defined s

\[\begin{split}\textbf{u}\cdot\nabla c = D\nabla^2c + R \\\end{split}\]

where \(c\) is the concentration field, \(\textbf{u}\) is the velocity, \(D\) is diffusivity, and \(R\) is the reaction.

Boundary conditions

Fixed value boundary condition

\[c = c_D \; \forall \textbf{x} \in \Gamma_D\]

Diffusive flux boundary condition

\[\nabla c \cdot \hat{n} = \textbf{h} \cdot \hat{n} \; \forall \textbf{x} \in \Gamma_N\]

where \(\hat{n}\) is the unit normal to \(\Gamma_N\)

Weak formulation

The weak formulation is stated as follows:

\[C := \{c \in H^1(\Omega) | c = c_D \;\text{on}\; \Gamma_D\}\]

and the corresponding test function

\[W := \{w \in H^1(\Omega) | w = 0 \;\text{on}\; \Gamma_D\}\]

The weak formulation, for \(c \in C\) and \(w \in W\)

\[0 = \left(w, \textbf{u}\cdot\nabla c\right)_\Omega + \left( \nabla w, D\nabla c \right)_\Omega - \left(w, R\right)_\Omega - \left(w, \textbf{h}\cdot \hat{n}\right)_{\Gamma_n} + \sum_{\Omega_e} \int_{\Omega_e} \tau \mathcal{R} \left(\textbf{u} \cdot \nabla w \right) d\Omega\]

where ther residue \(\mathcal{R}\) is

\[\mathcal{R} = \textbf{u} \cdot \nabla c - D\nabla^2c - R\]

where \(\mathcal{S}\) is the stabilization term for advection dominated problem. There are three types of \(\tau\) that the code provides. These are the shakib, su, codina type. For given a cell diameter \(h\), a cell Peclet number \(Pe_h=\frac{\lvert \textbf{u} \rvert h}{2D}\)

Stabilization parameters

Stabilization parameters \(\tau\) are predefined values via the get_stab_constant(tau_type) method. The tau_type parameter can be either shakib, su, or codina, and are defined as

\[\begin{split}\tau_{shakib} = \left( \left( \frac{2\lvert \textbf{u} \rvert^2}{h}\right)^2 + 9*\left(\frac{4D}{h}\right)^2 + R^2 \right)^{-0.5} \\\end{split}\]

\[\begin{split}\tau_{su} = \frac{h}{2 \lvert \textbf{u} \rvert} \left( coth(Pe_h) - \frac{1}{Pe_h}\right) \\\end{split}\]

\[\begin{split}\tau_{codina} = \frac{h}{2 \lvert \textbf{u} \rvert} \left(1 + \frac{1}{Pe_h} + \frac{hR}{2\lvert \textbf{u} \rvert}\right)^{-1} \\\end{split}\]

Transient Scalar-Transport

class flatiron_tk.physics.transient_scalar_transport.TransientScalarTransport(mesh, dt=0.01, theta=0.5, *args, **kwargs)

Bases: SteadyScalarTransport

add_stab()

Add stabilization term to the weak form.

build_function_space()

Build the function space for the transient scalar transport problem.

get_residual()

Compute the residual for the transient scalar transport problem.

Returns:

The residual expression for the transient scalar transport equations.

set_advection_velocity(u0, un)

Set the advection velocity for the transient scalar transport problem.

Parameters:

u0: The advection velocity at the previous time step. un: The advection velocity at the current time step.

set_diffusivity(D0, Dn)

Set the diffusivity for the transient scalar transport problem.

Parameters:

D0: The diffusivity at the previous time step. Dn: The diffusivity at the current time step.

set_initial_condition(u0)

Set the initial condition for the problem.

Parameters:

u0 (dolfinx.fem.Function) – The initial condition function.

set_midpoint_theta(theta)

Set the theta parameter for the implicit-explicit scheme.

Parameters:

theta: The theta parameter.

set_reaction(R0, Rn)

Set the reaction term for the transient scalar transport problem.

Parameters:

R0: The reaction term at the previous time step. Rn: The reaction term at the current time step.

set_time_step_size(dt)

Set the time step size for the transient scalar transport problem.

Parameters:

dt: The time step size.

set_weak_form(stab=False)

Set weak formulation as: dc/dt + mid(advection) - mid(diffusion) - mid(reaction) == 0

for a mid point type method where mid(f) = theta*f1 + (1-theta)*f0

Parameters:

stab (A boolean indicating whether to include stabilization terms in the weak form.)

stabilization_constant()

tau from Donea’s book Remark 5.8

update_previous_solution()

Update the previous solution with the current solution.

This class defines the transient version of the scalar transport problem

Strong formulation

The strong form of this now the time-dependent version of the TransientScalarTransport problem

\[\begin{split}\frac{\partial c}{\partial t} + \textbf{u}\cdot\nabla c = D\nabla^2c + R \\\end{split}\]

Initial and Boundary conditions

The boundary conditions are the same as the ones defined in the TransientScalarTransport class.

Here, the initial condition can be set through the set_initial_condition method.

Time discretization

In this module, the default time discretization is the midpoint \(\theta\) method defined as follows:

Let subscript \(()_0\) denote the variable in the previous time step and \(()_n\) define the variable at the current time step

Find \(c_n \in C\) and \(w \in W\) such that

\[0 = \left(w, \frac{c_n - c_0}{\Delta t} \right) + (1-\theta)\mathcal{L}(c_0; \textbf{u}_0, D_0, R_0) + \theta\mathcal{L}(c_n; \textbf{u}_n, D_n, R_n) + \sum_{\Omega_e} \int_{\Omega_e} \tau \mathcal{R} \left(\textbf{u} \cdot \nabla w \right) d\Omega\]

where

\[\mathcal{L}(c; \textbf{u}, D, R) := \left(w, \textbf{u}\cdot\nabla c\right)_\Omega + \left( \nabla w, D\nabla c \right)_\Omega - \left(w, R\right)_\Omega - \left(w, \textbf{h}\cdot \hat{n}\right)_{\Gamma_n}\]

where \(\mathcal{S}\) is the stabilization term for advection dominated problem. There are three types of \(\tau\) that the code provides. These are the shakib, su, codina type. For given a cell diameter \(h\), a cell Peclet number \(Pe_h=\frac{\lvert \textbf{u} \rvert h}{2D}\).

Stabilization parameters

Here, we only provide one stabilization constant for the SUPG stabilization based on the cordina stabilization parameter

\[\begin{split}\tau_{shakib} = \left( \left( \frac{1}{\theta \Delta t} \right)^2 + \left( \frac{2\lvert \textbf{u} \rvert^2}{h}\right)^2 + 9\left(\frac{4D}{h}\right)^2 + R^2 \right)^{-0.5} \\\end{split}\]

Class definition

class flatiron_tk.physics.TransientScalarTransport(mesh, dt=0.01, theta=0.5, *args, **kwargs)

add_stab()

Add stabilization term to the weak form.

build_function_space()

Build the function space for the transient scalar transport problem.

get_residual()

Compute the residual for the transient scalar transport problem.

Returns:

The residual expression for the transient scalar transport equations.

set_advection_velocity(u0, un)

Set the advection velocity for the transient scalar transport problem.

Parameters:

u0: The advection velocity at the previous time step. un: The advection velocity at the current time step.

set_diffusivity(D0, Dn)

Set the diffusivity for the transient scalar transport problem.

Parameters:

D0: The diffusivity at the previous time step. Dn: The diffusivity at the current time step.

set_initial_condition(u0)

Set the initial condition for the problem.

Parameters:

u0 (dolfinx.fem.Function) – The initial condition function.

set_midpoint_theta(theta)

Set the theta parameter for the implicit-explicit scheme.

Parameters:

theta: The theta parameter.

set_reaction(R0, Rn)

Set the reaction term for the transient scalar transport problem.

Parameters:

R0: The reaction term at the previous time step. Rn: The reaction term at the current time step.

set_time_step_size(dt)

Set the time step size for the transient scalar transport problem.

Parameters:

dt: The time step size.

set_weak_form(stab=False)

Set weak formulation as: dc/dt + mid(advection) - mid(diffusion) - mid(reaction) == 0

for a mid point type method where mid(f) = theta*f1 + (1-theta)*f0

Parameters:

stab (A boolean indicating whether to include stabilization terms in the weak form.)

stabilization_constant()

tau from Donea’s book Remark 5.8

update_previous_solution()

Update the previous solution with the current solution.