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()[source]¶

Print program name and information.

That is:

-> Version: 1.4
-> 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()[source]¶

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"

fenicsR13.geoToH5.geo_to_h5()[source]¶: 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)[source]¶

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
# - massflow: List of BC IDs to compute massflow J=int_bc dot(u,n) ds
# - 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
massflow: []
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)[source]¶: Construct the Input class.

get_from_input(map_list)[source]¶

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)[source]¶: 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)[source]¶

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)[source]¶

Construct the object.

This includes:

Read mesh

Read subdomains

Read boudnaries

fenicsR13.postprocessor module¶

Module to store the postprocessor classes.

class fenicsR13.postprocessor.Postprocessor(data, output_folder)[source]¶

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)[source]¶

Construct postprocessor.

Only sets arguments

plot_errors(show_popup)[source]¶

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()[source]¶

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)[source]¶

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)[source]¶: 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

__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()[source]¶

assemble()[source]¶

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()[source]¶

Solve the previously assembled system.

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"]))

Note

The following does not work in parallel because the mean is then only local. So convergence studies have to be performed in serial:

mean = np.mean(scalar_function.compute_vertex_values())

__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 –

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)[source]¶

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

Returns the same matrix if 2D

fenicsR13.tensoroperations.gen3DTFdim2(sigma)[source]¶

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)[source]¶

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)[source]¶

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)[source]¶

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)[source]¶

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)[source]¶

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)[source]¶

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)[source]¶

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)[source]¶

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}\]

fenicsR13 package¶

Submodules¶

fenicsR13.fenicsR13 module¶

fenicsR13.geoToH5 module¶

fenicsR13.input module¶

fenicsR13.meshes module¶

fenicsR13.postprocessor module¶

fenicsR13.solver module¶

fenicsR13.tensoroperations module¶

Module contents¶

fenicsR13

Navigation

Related Topics