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.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
orexamples
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:...
- 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
fenicsR13.meshes module¶
Module to store the mesh classes.
Currently, the only mesh format is h5.
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.
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
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
- 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 fileinput.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.
- 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
orr13
Heatflux
heat
orr13
Pressure
stress
orr13
Velocity
stress
orr13
Stress
stress
orr13
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 –
See also
calculate_errors
Function to return all field errors
Notes
For other norm types, see DOLFIN documentation 6 and search for norms.
References
- calculate_errors()[source]¶
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()[source]¶
Write all solver data to separate folder.
This includes the writing of:
Solutions
Parameter functions
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:
Heat source as f_mass
Mass Source as f_heat
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
- __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)[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}\]