Toggle navigation

.rst .pdf

Contents

fenicsR13 package

Submodules

fenicsR13.fenicsR13 module

Program to solve linearized R13 equations.

Different modes allow to solved the ecoupled heat or stress system. Different meshes can be used to perform a convergence study with given exact solution.

fenicsR13.fenicsR13.print_information()

Print program name and information.

That is:

-> Version: 1.5
-> Contact: Lambert Theisen <lambert.theisen@rwth-aachen.de>
-> Contact: Manuel Torrilhon <mt@mathcces.rwth-aachen.de>
-> Repository: <https://git.rwth-aachen.de/lamBOO/fenicsR13>
-> Documentation: <https://lamboo.pages.rwth-aachen.de/fenicsR13/>
  __            _          ____  _ _____
 / _| ___ _ __ (_) ___ ___|  _ \/ |___ /
| |_ / _ \ '_ \| |/ __/ __| |_) | | |_ \
|  _|  __/ | | | | (__\__ \  _ <| |___) |
|_|  \___|_| |_|_|\___|___/_| \_\_|____/

fenicsR13.fenicsR13.main()

Execute the main program.

Searches for an "input.yml" file in the current directory to use as input.

Usage:

# Install fenicsR13
pip install .

# Usage: <path_to_program> <input_file>
# Goto case folder:
cd tests/2d_r13
fenicsR13 inputs/           r13_1_coeffs_nosources_norot_inflow_p1p1p1p1p1_stab.yml

fenicsR13.geoToH5 module

Converter from geo-format to a mesh in h5-format.

Installation:

pip install .

Usage:

Usage: python3 <path_to_geoToH5.py> <geo_file> <h5_file>             [<gmsh cli arguments>]
E.g.: geoToH5 lid.geo lid5.h5 "-setnumber p 5"

3D Meshes: Append “-3” to the gmsh CLI arguments like:

geoToH5 mesh3d.geo mesh3d.h5 "-3"

fenicsR13.geoToH5.geo_to_h5()

Convert given geo-file to a h5-mesh.

fenicsR13.input module

Module for input related Classes.

Contains the Input class.

class fenicsR13.input.Input(yaml_file)

Bases: object

Class to handle the input file in YAML format.

An example input file could look like:

# General
# =======
# - output_folder: Used as output folder
output_folder: r13_1_coeffs_sources_rot_noinflow_p2p2p2p2p2_stab

# Meshes
# ======
# - meshes: List of input meshes in h5 format to run simulations on
meshes:
  - ../mesh/ring0.h5
  - ../mesh/ring1.h5
  - ../mesh/ring2.h5
  - ../mesh/ring3.h5
  - ../mesh/ring4.h5
  # - ../mesh/ring5.h5
  # - ../mesh/ring6.h5
  # - ../mesh/ring7.h5
  # - ../mesh/ring8.h5

# Numerical Parameters
# ====================
# - elements: Must contain the fields: theta, s, p, u, sigma
#   - fields: List of FEM parameters (shape, degree)
#     - shape: Element shape, e.g. Lagrange
#     - degree: Element degree, e.g. 2
# - stabilization: Must contain cip and gls
#   - cip: Collection of Continous Interior Penalty (CIP) parameters
#     - enable: Enable CIP stabilization
#     - delta_theta: Stabilization of grad(T)*grad(T_test) over edge
#     - delta_u: Stabilization of grad(u)*grad(u_test) over edge
#     - delta_p: Stabilization of grad(p)*grad(p_test) over edge
#   - gls: Collection of Garlerkin Least Squares (GLS) parameters
#     - enable: Enable GLS stabilization
#     - tau_energy: Stabilization with energy eq. residual
#     - tau_heatflux: Stabilization with heatflux eq. residual
#     - tau_mass: Stabilization with mass eq. residual
#     - tau_momentum: Stabilization with momentum eq. residual
#     - tau_stress: Stabilization with stress eq. residual
elements:
  theta:
    shape: Lagrange
    degree: 2
  s:
    shape: Lagrange
    degree: 2
  p:
    shape: Lagrange
    degree: 2
  u:
    shape: Lagrange
    degree: 2
  sigma:
    shape: Lagrange
    degree: 2
stabilization:
  cip:
    enable: True
    delta_theta: 1.0
    delta_u: 1.0
    delta_p: 0.01
  gls:
    enable: False
    tau_energy: 0.001
    tau_heatflux: 0.001
    tau_mass: 0.01
    tau_momentum: 0.01
    tau_stress: 0.01

# Formulation Parameters
# ======================
# - nsd: Number of spatial dimensions == 2
# - mode: Formulation mode, one of heat, stress, r13
# - kn: Knudsen numberkn
# - heat_source: Heat source function for mode==heat||r13
# - mass_source: Mass source function for mode==stress||r13
# - body_force: Body force for mode==stress||r13
nsd: 2
mode: r13
heat_source: 0
mass_source: 1.0 * (1.0 - (5.0*pow(R,2))/(18.0*pow(0.1,2))) * cos(phi)
body_force: [0,0]

# Region Parameters
# =================
# - regs: Dictionary of all mesh regions
#   - reg_id: Must contain the following parameters:
#     - kn: Knudsen number

# Boundary Conditions
# ===================
# - bcs: Dictionary of all boundary IDs from mesh
#   - bc_id: must contain the following parameters
#     - chi_tilde: Refaction coefficient in Maxwell accomodation model
#     - theta_w: Value for temperature at wall
#     - u_t_w: Value for tangential velocity at wall
#     - u_n_w: Value for normal velocity at wall
#     - p_w: Value for pressure at wall
#     - epsilon_w: Inflow-model parameter <=> Weight of pressure
bcs:
  3000:
    chi_tilde: 1.0
    theta_w: 1.0
    u_t_w: -10
    u_n_w: 0
    p_w: 0
    epsilon_w: 0
  3100:
    chi_tilde: 1.0
    theta_w: 0.5
    u_t_w: 0
    u_n_w: 0
    p_w: 0
    epsilon_w: 0

# Convergence Study
# =================
# - enable: Enable convergence study on given meshes
# - exact_solution: Path to exact solution in cpp-format to compare
# - plot: Show errors in matplotlib window. PDF output is default
# - write_systemmatrix: Writes out systemmatrix (LHS)
# - rescale_pressure: Rescale numerical pressure result for zero mean
# - relative_errors: Use relative errors. If esol=0, use absolute.
convergence_study:
  enable: True
  exact_solution: esols/1_coeffs_sources_rot_noinflow.cpp
  plot: False
  write_systemmatrix: False
  rescale_pressure: True
  relative_error: True

# Postprocessing
# ==============
# - write_pdfs: Write all solution fields as PDF plot
# - write_vecs: Write all solution fields as vectors
# - flows: List of BC IDs for <mass|heat> flow J=int_bc dot(<u|s>,n)dl
# - line_integrals: List of line integral dicts:
#   - name: Name for output
#   - expr: Expression to evaluate
#   - start: Start point
#   - end: End point
#   - res: Sampling resolution of line
postprocessing:
write_pdfs: True
write_vecs: True
flows: []
line_integrals:
  - name: "avg(abs(uy(x,y=0.5)))"
    expr: abs(uy)/8
    start: [0, 0.5]
    end: [8, 0.5]
    res: 10000
  - name: "avg(abs(uy(x,y=4.5)))"
    expr: abs(uy)/8
    start: [0, 4.5]
    end: [8, 4.5]
    res: 10000

# Parameter Study
# ==============
# - enable: Repeat simulation with different p. values (study)
# - parameter_key: Key as list, e.g. ["elemenets", "p", "degree"]
# - parameter_values: List of value for parameter, e.g. [0.01,0.1,1,10]
parameter_study:
enable: True
parameter_key: ["kn"]
parameter_values: [1,2,3]

Further input examples can be found in the tests or examples folders.

Examples

>>> corrupt_input = Input("etc/doctest_data/corrupt_input.yml")
Traceback (most recent call last):
...
Exception: Parsing error

>>> working_input = Input(
...     "tests/2d_heat/inputs/heat_01_coeffs_p1p1_stab.yml"
... )
Input:...

__init__(yaml_file)

Construct the Input class.

get_from_input(map_list)

Get value of the input dict by using a list of stacked keys.

Source: https://stackoverflow.com/questions/14692690/..

..access-nested-dictionary-items-via-a-list-of-keys/37704379

set_in_input(map_list, value)

Change value of the input dict by using a list of stacked keys.

fenicsR13.meshes module

Module to store the mesh classes.

Currently, the only mesh format is h5.

class fenicsR13.meshes.H5Mesh(h5_file)

Bases: object

Mesh class.

Raises:

Exception – File not found

Examples

>>> mesh = H5Mesh("non_existing_mesh.h5")
Traceback (most recent call last):
...
Exception: non_existing_mesh.h5 not found

__init__(h5_file)

Construct the object.

This includes:

1. Read mesh

2. Read subdomains

3. Read boudnaries

fenicsR13.postprocessor module

Module to store the postprocessor classes.

class fenicsR13.postprocessor.Postprocessor(data, output_folder)

Bases: object

The class for postprocessing.

The constructor expects a list of data dictionaries of the form:

[
    {
        'h': 0.9886573325052778,
        'theta': {'L_2': 0.18459230027019588,
        'l_inf': 0.11513941287387819},
        'sx': {'L_2': 0.38232086596749726,
        'l_inf': 0.30555017240112986},
        'sy': {'L_2': 0.51913814853123219,
        'l_inf': 0.2889203116681428}
    },
    {
        'h': 0.6340332990709842,
        'theta': {'L_2': 0.19952286856390053,
        'l_inf': 0.18948213107495288},
        'sx': {'L_2': 0.3528639452958619,
        'l_inf': 0.30118676712697767},
        'sy': {'L_2': 0.34778282823921497,
        'l_inf': 0.30792984640130788}
    }
]

Parameters:

• data (list of dicts) – Data to plot, see above for shape.

• output_folder (string) – output folder for PDF plot.

__init__(data, output_folder)

Construct postprocessor.

Only sets arguments

plot_errors(show_popup)

Use matplotlib to plot all errors in a figure.

A slope marker is added.

Exporting PDFs with:

# Write PDF without access to backend, i.e. use PDF backend
import matplotlib
matplotlib.use('pdf')
import matplotlib.pyplot as plt # pylint: disable=C0413

write_errors()

Write errors to a csv file.

The output file (e.g. error.csv) looks like:

h,p_L_2,p_l_inf,ux_L_2,ux_l_inf,uy_L_2,uy_l_inf,sigmaxx_L_2,...
0.9,0.2,0.2,0.1,0.1,0.2,0.1,0.2,0.5,0.2,0.3,0.2,0.5
0.6,0.1,0.1,0.0,0.1,0.1,0.1,0.0,0.3,0.0,0.3,0.0,0.3

fenicsR13.solver module

Solver module, contains the Solver class.

For usage examples, see the solver.Solver description.

class fenicsR13.solver.Solver(params, mesh, time)

Bases: object

Class to store the actual solver.

Possible order of methods in context of convergence study (see main program): “mesh=meshes[i=0]”, “__init__”, “assemble()”, “solve()”, “write()”, “…”, “mesh=meshes[i+1]”, “__init__”, “…”

Parameters:

• params (input.Input) – Input parameters

• mesh (meshes.H5Mesh) – Mesh

• time (string) – Identifier for e.g. convergence study naming

Returns:

Solver object

Return type:

Solver

Example

>>> # Example usage:
>>> from fenicsR13.input import Input
>>> from fenicsR13.meshes import H5Mesh
>>> params = Input(
...     "tests/2d_heat/inputs/heat_01_coeffs_p1p1_stab.yml"
... ) 
Input:...
>>> msh = H5Mesh("tests/2d_mesh/ring0.h5")
>>> solver = Solver(params.dict, msh, "0") # "0" means time=0

__init__(params, mesh, time)

Initialize solver and setup variables from input parameters.

params

Doctest

__createSolMacroScaExpr(cpp_string)

Return a DOLFIN scalar expression with predefined macros after solve.

These macros include:

Name	Macro CPP	Replacement
theta	theta	sol["theta"]
sx	sx	sol["s"].split()[0]
sy	sy	sol["s"].split()[1]
p	p	sol["p"]
ux	ux	sol["u"].split()[0]
uy	uy	sol["u"].split()[1]
sigmaxx	sigmaxx	sol["sigma"].split()[0]
sigmaxy	sigmaxy	sol["sigma"].split()[1]
sigmayx	sigmayx	sol["sigma"].split()[1]
sigmayy	sigmayy	sol["sigma"].split()[2]

__createMacroScaExpr(cpp_string)

Return a DOLFIN scalar expression with predefined macros.

These macros include:

Name	Macro	CPP Replacement
Radius wrt. to \((0,0)\)	R	sqrt(pow(x[0],2)+pow(x[1],2))
Angle wrt. \((0,0)\)	phi	atan2(x[1],x[0])

The following expressions are therefore equal:

# expr1 is equal to expr2
expr1 = self.__createMacroScaExpr("R*cos(phi)")
expr2 = dolfin.Expression(
    "R*cos(phi)",
    degree=2,
    R=dolfin.Expression("sqrt(pow(x[0],2)+pow(x[1],2))", degree=2),
    phi=dolfin.Expression("atan2(x[1],x[0])", degree=2),
)

__createMacroVecExpr(cpp_strings)

Return a DOLFIN vector expression with predefined macros.

These macros include:

Name	Macro	CPP Replacement
Radius wrt. to \((0,0)\)	R	sqrt(pow(x[0],2)+pow(x[1],2))
Angle wrt. \((0,0)\)	phi	atan2(x[1],x[0])

See also

_Solver__createMacroScaExpr

__createMacroTenExpr(cpp_strings)

Return a DOLFIN tensor expression with predefined macros.

These macros include:

Name	Macro	CPP Replacement
Radius wrt. to \((0,0)\)	R	sqrt(pow(x[0],2)+pow(x[1],2))
Angle wrt. \((0,0)\)	phi	atan2(x[1],x[0])

See also

_Solver__createMacroScaExpr

__setup_function_spaces()

Set up function spaces for trial and test functions for assembling.

Depends on the mode. Function spaces depend on the choice of the element and its degree (see the input file input.Input).

The following DOLFIN functions are used:

field	DOLFIN Function
theta	FiniteElement
s	VectorElement
p	FiniteElement
u	VectorElement
sigma	TensorElement

__check_regions()

Check if all regions from the input mesh have params prescribed.

Raises an exception if one region is missing.

__check_bcs()

Check if all boundaries from the input mesh have BCs prescribed.

Raises an exception if one BC is missing.

normal()

assemble()

Assemble the weak form of the system, depending on the mode. This routine contains the main implementation of a stable variational formulation of the R13 system which allows equal order ansatz spaces for all variables and correctly incorporates the R13 boundary conditions.

The system of partial differential equations results from the two dimensional, linearized R13 equations of [1], that have been derived from Boltzmann’s equation as a continuum model for rarefied gas flows. The boundary conditions have been presented in [2], while the variational formulation and its stabilization was introduced in [3]. See also [4] for a general review and [5] for the modified stable boundary conditions.

[1]

H. Struchtrup, M. Torrilhon (2003). Regularization of Grad’s 13 Moment Equations: Derivation and Linear Analysis.

[2]

M. Torrilhon, H. Struchtrup (2007). Boundary Conditions for Regularized 13-Moment-Equations for Micro-Channel-Flows

[3]

A. Westerkamp, M. Torrilhon (2019). Finite Element Methods for the Linear Regularized 13-Moment Equations Describing Slow Rarefied Gas Flows

[4]

M. Torrilhon (2016). Modeling Nonequilibrium Gas Flow Based on Moment Equations

[5]

M. Torrilhon, N. Sarna (2017). Hierarchical Boltzmann Simulations and Model Error Estimation

Identities:

• \(\boldsymbol{x}_1 \in \mathbb{R}^3\), \(\boldsymbol{x}_2 \in \mathbb{R}^{3 \times 3}\), …

• Test function \(\boldsymbol{\psi}\) is assumed to be symmetric and trace-less:

\[\begin{split}\boldsymbol{\psi} : \boldsymbol{I} = 0 \\ (\boldsymbol{\psi})_{\text{sym}} = \boldsymbol{\psi}\end{split}\]

• Trace of vector gradient is divergence of vector:

\[\langle \boldsymbol{I},\boldsymbol\nabla \boldsymbol{x}_1 \rangle = \boldsymbol{I} : \boldsymbol\nabla \boldsymbol{x}_1 = \text{div}(\boldsymbol{x}_1)\]

• Inner product has orthogonality property with respect to the additive symmetric/skewsymmetric tensor decomposition:

\[\begin{split}\langle (\boldsymbol{x}_2)_{\text{sym}} , \boldsymbol{y}_2 \rangle &= \langle (\boldsymbol{x}_2)_{\text{sym}} , (\boldsymbol{y}_2)_{\text{sym}} + (\boldsymbol{y}_2)_{\text{skew}} \rangle \\ &= \langle (\boldsymbol{x}_2)_{\text{sym}} , (\boldsymbol{y}_2)_{\text{sym}} \rangle + \langle (\boldsymbol{x}_2)_{\text{sym}} , (\boldsymbol{y}_2)_{\text{skew}} \rangle \\ &= \langle (\boldsymbol{x}_2)_{\text{sym}} , (\boldsymbol{y}_2)_{\text{sym}} \rangle\end{split}\]

• Inner product of STF tensor with arbitrary tensor:

\[T_{i_{1} i_{2} \ldots i_{n} k_{1} \cdots k_{r}} A_{i_{1} i_{2} \ldots i_{n} j_{1} \cdots j_{n}} = T_{i_{1} i_{2} \ldots i_{n} k_{1} \cdots k_{r}} A_{\langle i_{1} i_{2} \ldots i_{n}\rangle j_{1} \cdots j_{n}}\]

Tricks of the trade:

\[\begin{split}\langle {(\boldsymbol\nabla \boldsymbol{s})}_{\text{STF}} , \boldsymbol\nabla \boldsymbol{r} \rangle &= \langle {(\boldsymbol\nabla \boldsymbol{s})}_{\text{sym}} - \frac13\text{tr}(\boldsymbol\nabla \boldsymbol{s})\boldsymbol{I} , \boldsymbol\nabla \boldsymbol{r} \rangle \\ &= \langle {(\boldsymbol\nabla \boldsymbol{s})}_{\text{sym}} , \boldsymbol\nabla \boldsymbol{r} \rangle -\frac13\text{tr}(\boldsymbol\nabla \boldsymbol{s}) \langle \boldsymbol{I},(\boldsymbol\nabla \boldsymbol{r}) \rangle \\ &= \langle {(\boldsymbol\nabla \boldsymbol{s})}_{\text{sym}} , {(\boldsymbol\nabla \boldsymbol{r})}_{\text{sym}} \rangle -\frac13\text{tr}(\boldsymbol\nabla \boldsymbol{s}) \text{tr}(\boldsymbol\nabla \boldsymbol{r}) \\ &= \langle {(\boldsymbol\nabla \boldsymbol{s})}_{\text{sym}} , {(\boldsymbol\nabla \boldsymbol{r})}_{\text{sym}} \rangle -\frac13\text{div}(\boldsymbol{s}) \text{div}(\boldsymbol{r})\end{split}\]

\[\begin{split}\langle {(\boldsymbol\nabla \boldsymbol{s})}_{\text{STF}} , \boldsymbol{\psi} \rangle &= \langle {(\boldsymbol\nabla \boldsymbol{s})}_{\text{sym}} - \frac13\text{tr}(\boldsymbol\nabla \boldsymbol{s})\boldsymbol{I} , \boldsymbol{\psi} \rangle \\ &= \langle {(\boldsymbol\nabla \boldsymbol{s})}_{\text{sym}} , \boldsymbol{\psi} \rangle -\frac13\text{tr}(\boldsymbol\nabla \boldsymbol{s}) \langle \boldsymbol{I},\boldsymbol{\psi} \rangle \\ &= \langle {\boldsymbol\nabla \boldsymbol{s}} , {\boldsymbol{\psi}} \rangle -\frac13\text{div}(\boldsymbol{s}) \text{tr}(\boldsymbol{\psi}) \\ &= \langle {\boldsymbol\nabla \boldsymbol{s}} , {\boldsymbol{\psi}} \rangle\end{split}\]

solve()

Solve the previously assembled system.

Note that setting the mg levels for gamg doesn’t work (see here_). Also see some note here2_.

_here: https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2019-April/038223.html # noqa _here2: https://www.firedrakeproject.org/demos/saddle_point_systems.py.html # noqa

Some available solver options:

# Some solver params
solver_parameters={
    "linear_solver": "gmres", "preconditioner": "ilu" # or
    "linear_solver": "petsc", "preconditioner": "ilu" # or
    "linear_solver": "direct" # or
    "linear_solver": "mumps" # or
    "linear_solver": "mumps"
}

# List all available solvers and preconditioners:

list_linear_solver_methods()
bicgstab  | Biconjugate gradient stabilized method
cg        | Conjugate gradient method
default   | default linear solver
gmres     | Generalized minimal residual method
minres    | Minimal residual method
mumps     | MUMPS (MUltifrontal Massively Parallel Sparse Direct)
petsc     | PETSc built in LU solver
richardson| Richardson method
superlu   | SuperLU
tfqmr     | Transpose-free quasi-minimal residual method
umfpack   | UMFPACK (Unsymmetric MultiFrontal sparse LU factoriz.)
# "direct" means "default" means "lu" of default backend
print(parameters["linear_algebra_backend"]) # usually PETSc

list_krylov_solver_preconditioners()
amg              |  Algebraic multigrid
default          |  default preconditioner
hypre_amg        |  Hypre algebraic multigrid (BoomerAMG)
hypre_euclid     |  Hypre parallel incomplete LU factorization
hypre_parasails  |  Hypre parallel sparse approximate inverse
icc              |  Incomplete Cholesky factorization
ilu              |  Incomplete LU factorization
jacobi           |  Jacobi iteration
none             |  No preconditioner
petsc_amg        |  PETSc algebraic multigrid
sor              |  Successive over-relaxation

__load_exact_solution()

Load exact solution from the location given in input.yml.

The exact solution must be C++ format with a specific syntax. The esol.cpp must contain the classes:

Class	mode
Temperature	heat or r13
Heatflux	heat or r13
Pressure	stress or r13
Velocity	stress or r13
Stress	stress or r13

The file has to follow a specific syntax for DOLFIN. An example file could look like:

#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
#include <cmath>
// additional includes
#include <boost/math/special_functions/bessel.hpp>
using namespace std;
namespace py = pybind11;
#include <dolfin/function/Expression.h>
double lambda_3 = sqrt(3.0/2.0); // some own constants
class Temperature : public dolfin::Expression {
    public:
    Temperature() : dolfin::Expression() {}
    void eval(Eigen::Ref<Eigen::VectorXd> values,
            Eigen::Ref<const Eigen::VectorXd> x) const override {
        values[0] = 1; // value
    }
};
class Heatflux : public dolfin::Expression {
    public:
    Heatflux() : dolfin::Expression(2) {} // note components=2!
    void eval(Eigen::Ref<Eigen::VectorXd> values,
            Eigen::Ref<const Eigen::VectorXd> x) const override {
        values[0] = 42;
        values[1] = 3.14;
    }
};
class Pressure : public dolfin::Expression {
    public:
    Pressure() : dolfin::Expression() {}
    void eval(Eigen::Ref<Eigen::VectorXd> values,
            Eigen::Ref<const Eigen::VectorXd> x) const override {
        values[0] = boost::math::cyl_bessel_i(1,2.71); // external
    }
};
class Velocity : public dolfin::Expression {
    public:
    Velocity() : dolfin::Expression(2) {}
    void eval(Eigen::Ref<Eigen::VectorXd> values,
            Eigen::Ref<const Eigen::VectorXd> x) const override {
        values[0] = lambda_3;
        values[1] = 2;
    }
};
class Stress : public dolfin::Expression {
    public:
    Stress() : dolfin::Expression(2,2) {} // note dim=2, shape=(2,2)
    void eval(Eigen::Ref<Eigen::VectorXd> values,
            Eigen::Ref<const Eigen::VectorXd> x) const override {
        double xx_val = 1.23;
        double xy_val = 1.23;
        double yy_val = 1.23;
        values[0] = xx_val;
        values[1] = xy_val;
        values[2] = yy_val;
        // values[3] = xy_val // not used due to symmetry, skip
    }
};
PYBIND11_MODULE(SIGNATURE, m) { // needed for DOLFIN
    py::class_<Temperature, std::shared_ptr<Temperature>,
               dolfin::Expression>
        (m, "Temperature")
    .def(py::init<>());
    py::class_<Heatflux, std::shared_ptr<Heatflux>,
               dolfin::Expression>
        (m, "Heatflux")
    .def(py::init<>());
    py::class_<Pressure, std::shared_ptr<Pressure>,
               dolfin::Expression>
        (m, "Pressure")
    .def(py::init<>());
    py::class_<Velocity, std::shared_ptr<Velocity>,
               dolfin::Expression>
        (m, "Velocity")
    .def(py::init<>());
    py::class_<Stress, std::shared_ptr<Stress>, dolfin::Expression>
        (m, "Stress")
    .def(py::init<>());
}

__calc_sf_mean(scalar_function)

Calculate the mean of a scalar function.

np.set_printoptions(precision=16)
# Precision is not soo nice, only 9 digits:
print(mean)
# In solve() has m. prec. hmmm:
print(self.__calc_sf_mean(self.sol["p"]))

__calc_sf_min(scalar_function)

Calculate the minimum of a scalar function (works in parallel).

__calc_field_errors(field_, field_e_, v_field, name_)

Calculate both \(L_2\) and \(l_\infty\) errors.

Works for scalars, vectors and tensors. The difference is written to a file. The exact solution is written to a file. Relative errors are based per field component and the maximum value of the analytical solution. If the analytical solution is uniformly zero, then the absolute erorrs is used. (equivalent to setting the maximum to 1)

Parameters:

• field (DOLFIN function) – calculated field

• field_e (DOLFIN function) – exact solution of field

• v_field (DOLFIN fspace) – function space for error calculation

• name (string) – name of the field, used to write difference

Returns:

Dict with an error list for “L_2” and a list for “l_inf”

Return type:

dict

Raises:

Nothing –

See also

calculate_errors

Function to return all field errors

Notes

For other norm types, see DOLFIN documentation [6] and search for norms.

The following is very slow, I don’t know why: .. code-block:: python

err_f_L2_full = df.errornorm(

field_e_i, field_i, “L2”, degree_rise=1

)

References

[6]

DOLFIN documentation

calculate_errors()

Calculate and return the errors of numerical to exact solution.

This includes all calculated fields.

Note

Usage of np.max() does not work in parallel. So convergence studies have to be performed in serial for now. Final fields should be right, so MPI can be used for production simulations.

Returns:

dict – Errors

write_content_to_file(filename, content)

Write content to a file in the output folder.

write()

Write all solver data to separate folder.

This includes the writing of:

1. Solutions

2. Parameter functions

3. System matrices if set in input file

__write_solutions()

Write all solutions fields.

__write_vecs()

Write all solutions to vectors.

__write_parameters()

Write parameter functions for debug reasons.

This includes:

1. Heat source as f_mass

2. Mass Source as f_heat

3. Body force as f_body

The parameter functions are internally interpolated into a \(P_1\) space before writing.

__write_discrete_system()

Write the discrete system matrix and the RHS vector.

Can be used to analyze e.g. condition number. Include writing of \(\mathbf{A}\) and \(\mathbf{b}\) of \(\mathbf{A} \mathbf{x} = \mathbf{b}\).

File-ending is .mat.

Import the matrices/vectors e.g. into Matlab with:

% Input into MATLAB
A = table2array(readtable("A_0.mat","FileType","text"));
b = table2array(readtable("b_0.mat","FileType","text"));

Julia:

using DelimitedFiles
A = readdlm("A_0.mat", ' ', Float64, '\n')

Example

>>> # Construct LHS
>>> from dolfin import *
>>> mesh = IntervalMesh(2 ,0, 1)
>>> V = FunctionSpace(mesh, "Lagrange", 1)
>>> u = TrialFunction(V)
>>> v = TestFunction(V)
>>> a = inner(grad(u),grad(v))*dx
>>> L = df.Constant(1)*v*dx
>>> lhs = assemble(a)
>>> print(lhs.array())
[[ 2. -2.  0.]
 [-2.  4. -2.]
 [ 0. -2.  2.]]
>>> rhs = assemble(L)
>>> print(rhs.get_local())
[ 0.25  0.5   0.25]
>>> # Assign LHS to solver
>>> from fenicsR13.input import Input
>>> from fenicsR13.meshes import H5Mesh
>>> params = Input(
...     "tests/2d_heat/inputs/heat_01_coeffs_p1p1_stab.yml"
... ) 
Input:...
>>> msh = H5Mesh("tests/2d_mesh/ring0.h5")
>>> solver = Solver(params.dict, msh, "0") # "0" means time=0
>>> solver.form_lhs = a
>>> solver.form_rhs = L
>>> solver.output_folder = "./"
>>> solver._Solver__write_discrete_system()
Write ./A_0.mat
Write ./b_0.mat
>>> print(open("A_0.mat","r").read())
2.0000...00000e+00 -2.0000...00000e+00 0.000...000000e+00
-2.0000...00000e+00 4.0000...00000e+00 -2.0000...00000e+00
0.0000...00000e+00 -2.0000...00000e+00 2.000...000000e+00

>>> print(open("b_0.mat","r").read())
2.500000000000000000e-01
5.000000000000000000e-01
2.500000000000000000e-01

__write_xdmf(name, field, write_pdf)

Write a given field to a XDMF file in the output folder.

Arguments

name

The name to be given to the field. Will be used as filename and is the name of the field in e.g. Paraview.

field

The field to write.

write_pdf

If true, write a simple PDF plot for all solution fields

Problems

For MPi execution, XDMFFile.write function fails if the mesh is too small. I wasn’t able to fix this in 2 days. Possible reasons: - Too much overlap in Schwarz-like decomposition

__write_vec(name, field)

Write a given field to a vector.

fenicsR13.tensoroperations module

Module to gather all additional tensor operations not present in UFL.

This especially includes all 3D operations and operations on tensors with rank above 2.

[STR2005] (1,2)

Struchtrup, H. (2005). Macroscopic transport equations for rarefied gas flows. In Macroscopic Transport Equations for Rarefied Gas Flows. Springer, Berlin, Heidelberg.

[TOR2018] (1,2)

Torrilhon, M. et al. (2018). “Kinetic Theory of Non-Equilibrium Gas Flows: Theory and Computations”. Lecture Notes. Indian Institute of Technology Madras.

fenicsR13.tensoroperations.gen3DTFdim3(sigma)

Return the symmetric and trace-free matrix for a 2-rank tensor if 3D.

Returns the same matrix if 2D

fenicsR13.tensoroperations.gen3DTFdim2(sigma)

Return the symmetric and trace-free matrix for a 2-rank tensor if 2D.

Returns the same matrix if 3D

fenicsR13.tensoroperations.stf3d2(rank2_2d)

Return the 3D symmetric and trace-free part of a 2D 2-tensor.

Warning

Return a 2-tensor with the same dimensions as the input tensor.

For the \(2 \times 2\) case, return the 3D symmetric and trace-free (dev(sym(.))) \(B \in \mathbb{R}^{2 \times 2}\) of the 2D 2-tensor \(A \in \mathbb{R}^{2 \times 2}\).

\[\begin{split}A &= \begin{pmatrix} a_{xx} & a_{xy} \\ a_{yx} & a_{yy} \end{pmatrix} \\ B &= (A)_\mathrm{dev} = \frac{1}{2} (A)_\mathrm{sym} - \frac{1}{3} \mathrm{tr}(A) I_{2 \times 2}\end{split}\]

fenicsR13.tensoroperations.sym3d3(rank3_3d)

Return the symmetric part of a 3D 3-tensor.

Returns the symmetric part \(A_{(i j k)}\) of a 3D rank-3 tensor \(A \in \mathbb{R}^{3 \times 3 \times 3}\) [STR2005] [TOR2018].

\[A_{(i j k)}=\frac{1}{6}\left( A_{i j k}+A_{i k j}+A_{j i k}+A_{j k i}+A_{k i j}+A_{k j i} \right)\]

fenicsR13.tensoroperations.stf3d3(rank3_3d)

Return the symmetric and trace-free part of a 3D 3-tensor.

Return the symmetric and trace-free part of a 3D 3-tensor. (dev(sym(.)))

Returns the STF part \(A_{\langle i j k\rangle}\) of a rank-3 tensor \(A \in \mathbb{R}^{3 \times 3 \times 3}\) [TOR2018].

\[A_{\langle i j k\rangle}=A_{(i j k)}-\frac{1}{5}\left[A_{(i l l)} \delta_{j k}+A_{(l j l)} \delta_{i k}+A_{(l l k)} \delta_{i j}\right]\]

A gradient \(\frac{\partial S_{\langle i j}}{\partial x_{k \rangle}}\) of a symmetric 2-tensor \(S \in \mathbb{R}^{3 \times 3}\) has the deviator [STR2005]:

\[\begin{split}\frac{\partial S_{\langle i j}}{\partial x_{k \rangle}} =& \frac{1}{3}\left( \frac{\partial S_{i j}}{\partial x_{k}} +\frac{\partial S_{i k}}{\partial x_{j}} +\frac{\partial S_{j k}}{\partial x_{i}} \right) -\frac{1}{15} \biggl[ \left( \frac{\partial S_{i r}}{\partial x_{r}} +\frac{\partial S_{i r}}{\partial x_{r}} +\frac{\partial S_{r r}}{\partial x_{i}} \right) \delta_{j k} \\ &+ \left( \frac{\partial S_{j r}}{\partial x_{r}} +\frac{\partial S_{j r}}{\partial x_{r}} +\frac{\partial S_{r r}}{\partial x_{j}} \right) \delta_{i k} +\left( \frac{\partial S_{k r}}{\partial x_{r}} +\frac{\partial S_{k r}}{\partial x_{r}} +\frac{\partial S_{r r}}{\partial x_{k}} \right) \delta_{i j} \biggr]\end{split}\]

fenicsR13.tensoroperations.div3d3(rank3_3d)

Return the 3D divergence of a 3-tensor.

Return the 3D divergence of a 3-tensor as

\[{(\mathrm{div}(m))}_{ij} = \frac{\partial m_{ijk}}{\partial x_k}\]

Compare with [SCH2009] (p. 92).

[SCH2009]

Heinz Schade, Klaus Neemann (2009). “Tensoranalysis”. 2. überarbeitete Auflage.

fenicsR13.tensoroperations.gen3dTF2(rank2_2d)

Generate a 3D tracefree 2-tensor from a 2D 2-tensor.

Returns the synthetic 3D version \(A \in \mathbb{R}^{3 \times 3}\) of a 2D rank-2 tensor \(B \in \mathbb{R}^{2 \times 2}\).

\[\begin{split}B &= \begin{pmatrix} b_{xx} & b_{xy} \\ b_{yx} & b_{yy} \end{pmatrix} \\ A &= \begin{pmatrix} b_{xx} & b_{xy} & 0 \\ b_{yx} & b_{yy} & 0 \\ 0 & 0 & -(b_{yx}+b_{yy}) \end{pmatrix}\end{split}\]

Example

>>> t2d = df.Expression((("1","2"),("3","4")), degree=1)
>>> t3d = gen3dTF2(t2d)
>>> print(
...     t3d[0,0]==t2d[0,0] and
...     t3d[0,1]==t2d[0,1] and
...     t3d[0,2]==0 and
...     t3d[1,0]==t2d[1,0] and
...     t3d[1,1]==t2d[1,1] and
...     t3d[1,2]==0 and
...     t3d[2,0]==0 and
...     t3d[2,1]==0 and
...     t3d[2,2]==-t2d[0,0]-t2d[1,1]
... )
True

fenicsR13.tensoroperations.gen3d1(rank1_2d)

Return synthetic 3d version of 2d vector with zero last component.

Return synthetic 3d version of 2d vector \(s_{\mathrm{in}} = (s_x, s_y)\) as \(s_{\mathrm{out}} = (s_x, s_y, 0)\) .

fenicsR13.tensoroperations.gen3d2(rank2_2d)

Generate a 3D 2-tensor from a 2D 2-tensor (add zeros to last dimensions).

Returns the synthetic 3D version \(A \in \mathbb{R}^{3 \times 3}\) of a 2D rank-2 tensor \(B \in \mathbb{R}^{2 \times 2}\). The 3D-components are set to zero.

\[\begin{split}B &= \begin{pmatrix} b_{xx} & b_{xy} \\ b_{yx} & b_{yy} \end{pmatrix} \\ A &= \begin{pmatrix} b_{xx} & b_{xy} & 0 \\ b_{yx} & b_{yy} & 0 \\ 0 & 0 & 0 \end{pmatrix}\end{split}\]

Example

>>> t2d = df.Expression((("1","2"),("3","4")), degree=1)
>>> t3d = gen3d2(t2d)
>>> print(
...     t3d[0,0]==t2d[0,0] and
...     t3d[0,1]==t2d[0,1] and
...     t3d[0,2]==0 and
...     t3d[1,0]==t2d[1,0] and
...     t3d[1,1]==t2d[1,1] and
...     t3d[1,2]==0 and
...     t3d[2,0]==0 and
...     t3d[2,1]==0 and
...     t3d[2,2]==0
... )
True

fenicsR13.tensoroperations.grad3dOf2(rank2_3d, dim)

Return 3D gradient of 3D 2-tensor.

Returns the 3D gradient (w.r.t. \(x,y,z\)) \(\frac{\partial A_{ij}}{\partial x_{k}}\) of a 3D \(z\)-independent rank-2 tensor \(A(x,y) \in \mathbb{R}^{3 \times 3}\) where

\[\begin{split}A_{ij} = \begin{pmatrix} a_{xx}(x,y) & a_{xy}(x,y) & a_{xz}(x,y) \\ a_{yx}(x,y) & a_{yy}(x,y) & a_{yz}(x,y) \\ a_{zx}(x,y) & a_{zy}(x,y) & a_{zz}(x,y) \end{pmatrix}\end{split}\]

The function then returns the 3-tensor with components \(\frac{\partial A_{ij}}{\partial x_{k}}\) where

\[\begin{split}\frac{\partial A_{ij}}{\partial x_{1}} = \frac{\partial A_{ij}}{\partial x} &= \begin{pmatrix} \partial_x a_{xx}(x,y) & \partial_x a_{xy}(x,y) & \partial_x a_{xz}(x,y) \\ \partial_x a_{yx}(x,y) & \partial_x a_{yy}(x,y) & \partial_x a_{yz}(x,y) \\ \partial_x a_{zx}(x,y) & \partial_x a_{zy}(x,y) & \partial_x a_{zz}(x,y) \end{pmatrix} , \\ \frac{\partial A_{ij}}{\partial x_{2}} = \frac{\partial A_{ij}}{\partial y} &= \begin{pmatrix} \partial_y a_{xx}(x,y) & \partial_y a_{xy}(x,y) & \partial_y a_{xz}(x,y) \\ \partial_y a_{yx}(x,y) & \partial_y a_{yy}(x,y) & \partial_y a_{yz}(x,y) \\ \partial_y a_{zx}(x,y) & \partial_y a_{zy}(x,y) & \partial_y a_{zz}(x,y) \end{pmatrix} , \\ \frac{\partial A_{ij}}{\partial x_{3}} = 0 &= \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\end{split}\]

Module contents

previous

fenicsR13

next

2d_heat