From 0f6df1c248c693e73ae5c029cea11240c3800606 Mon Sep 17 00:00:00 2001 From: Ouardghi Date: Mon, 18 May 2026 16:47:55 +0200 Subject: [PATCH 1/2] =?UTF-8?q?Add=20Taylor=E2=80=93Green=20NS=20benchmark?= =?UTF-8?q?=20studying=20order=20reduction=20induced=20by=20time-dependent?= =?UTF-8?q?=20bcs?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../StroemungsRaum_FEniCSx/README.rst | 57 ++ .../StroemungsRaum_FEniCSx/environment.yml | 23 + ...rStokes_2D_TaylGreen_monolithic_FEniCSx.py | 668 ++++++++++++++++++ ..._equations_TaylGreen_monolithic_FEniCSx.py | 169 +++++ ...ier_Stokes_TaylGreen_monolithic_FEniCSx.py | 154 ++++ 5 files changed, 1071 insertions(+) create mode 100644 pySDC/projects/StroemungsRaum_FEniCSx/README.rst create mode 100644 pySDC/projects/StroemungsRaum_FEniCSx/environment.yml create mode 100755 pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py create mode 100755 pySDC/projects/StroemungsRaum_FEniCSx/run_Navier_Stokes_equations_TaylGreen_monolithic_FEniCSx.py create mode 100644 pySDC/projects/StroemungsRaum_FEniCSx/tests/test_Navier_Stokes_TaylGreen_monolithic_FEniCSx.py diff --git a/pySDC/projects/StroemungsRaum_FEniCSx/README.rst b/pySDC/projects/StroemungsRaum_FEniCSx/README.rst new file mode 100644 index 0000000000..092a077309 --- /dev/null +++ b/pySDC/projects/StroemungsRaum_FEniCSx/README.rst @@ -0,0 +1,57 @@ +StroemungsRaum_FEniCSx +====================== + +**StroemungsRaum_FEniCSx** is a research software project derived from the +original *StroemungsRaum* framework and developed within the context of the +BMBF-funded project + +*“StrömungsRaum – Novel Exascale Architectures with Heterogeneous Hardware +Components for Computational Fluid Dynamics Simulations”* +(October 2022 – September 2025). + +This repository provides a modernized implementation of the original +StroemungsRaum research codes using **FEniCSx** instead of the legacy FEniCS +framework. The migration enables compatibility with current finite element +software ecosystems and improved support for modern HPC architectures. + +Scope of This Repository +------------------------ +This repository contains the **Forschungszentrum Jülich (FZJ)** contribution to +the StrömungsRaum project implemented with FEniCSx, focusing on: + +- Parallel-in-time methods +- Combined space–time parallelization +- Scalable solvers for time-dependent PDEs +- High-performance finite element implementations with PETSc and MPI +- Research on exascale-ready CFD algorithms + +The primary objective is to expose concurrency beyond spatial parallelism and +enable efficient large-scale simulations on modern heterogeneous HPC systems. + +Software Stack +-------------- +This project is based on the modern FEniCSx ecosystem, including: + +- FEniCSx +- PETSc +- mpi4py +- petsc4py +- UFL +- Basix +- DOLFINx + +Research Focus +-------------- +The software serves as a research platform for investigating: + +- Scalable time integration methods +- Space–time parallelism +- High-order temporal discretizations +- Solver robustness for incompressible flows +- Parallel performance on heterogeneous architectures + +Funding +------- +Funded by the **German Federal Ministry of Education and Research (BMBF)** under +grant number **16ME0708**. + diff --git a/pySDC/projects/StroemungsRaum_FEniCSx/environment.yml b/pySDC/projects/StroemungsRaum_FEniCSx/environment.yml new file mode 100644 index 0000000000..45a1306206 --- /dev/null +++ b/pySDC/projects/StroemungsRaum_FEniCSx/environment.yml @@ -0,0 +1,23 @@ +--- + +name: pySDC +channels: + - conda-forge +dependencies: + + - fenics-dolfinx=0.10.0 + - ufl=2025.2.1 + - basix=0.10.0 + - mpi4py + - petsc4py + - slepc4py + - pyvista + - numpy + - scipy>=0.17.1 + - dill>=0.2.6 + - matplotlib>=3.0 + - pytest + - pytest-cov + - pip + - pip: + - qmat>=0.1.8 diff --git a/pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py b/pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py new file mode 100755 index 0000000000..9149e4bcb9 --- /dev/null +++ b/pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py @@ -0,0 +1,668 @@ +import numpy as np +from mpi4py import MPI +from petsc4py import PETSc +import ufl +import basix +import dolfinx as dfx +import dolfinx.fem.petsc +import dolfinx_mpc + +from pySDC.core.problem import Problem +from pySDC.implementations.datatype_classes.mesh import mesh + + +class fenicsx_NSE_mass(Problem): + r""" + Monolithic incompressible Navier–Stokes problem in 2D using FEniCSx. + + This class implements a finite element discretization of the incompressible + Navier–Stokes equations in a square domain using a monolithic formulation + for velocity and pressure. + + Governing equations: + .. math:: + \frac{\partial \mathbf{u}}{\partial t} + + (\mathbf{u} \cdot \nabla)\mathbf{u} + - \nu \Delta \mathbf{u} + + \nabla p = f, + \quad \nabla \cdot \mathbf{u} = 0 + + Domain: + A square domain :math:`\Omega = [-0.5, 0.5]^2` discretized with triangular elements. + + Spatial discretization: + - Continuous Lagrange elements for velocity and pressure + - Mixed finite element space :math:`W = V \times Q` + - Monolithic formulation for the coupled velocity–pressure system + + Manufactured solution: + A smooth, analytical, time-dependent solution is given by: + ..math:: + u(x,y,t) = 1 - e^{-8\pi^{2}\nu t}\sin(2\pi(x-t))\sin(\pi y)\cos(\pi y) + v(x,y,t) = - e^{-8\pi^{2}\nu t} \cos\(2\pi(x-t)) \cos^{2}(\pi y) + p(x,y,t) = 1 + \frac{4}{17} e^{-16\pi^{2}\nu t} \cos\(4\pi(x-t)) \cos(\pi y) + + The solution is specifically constructed to investigate the impact of boundary + conditions on numerical accuracy. The source term g =(g_1, g_2) is derived from + the exact solution to ensure that it satisfies the Navier–Stokes equations. + + Boundary conditions: + - Time-dependent Dirichlet conditions on the left and right boundaries + - Time-independent Dirichlet conditions on the top and bottom boundaries + + This combination is intentionally chosen to study order reduction effects + induced by time-dependent boundary conditions. + + In the related class fenicsx_NSE_periodic_mass, the same problem is formulated with + periodic boundary conditions on the left and right boundaries, while retaining the + time-independent Dirichlet conditions on the top and bottom boundaries, in order to + isolate and compare the effects of time-dependent boundary data on the numerical accuracy. + + + Parameters + ---------- + nelems : int + Number of elements per spatial direction. + t0 : float + Initial time. + family : str + Finite element family (e.g., 'CG' for Lagrange). + order : int + Polynomial degree of the function spaces. + nu : float + Kinematic viscosity. + comm : MPI.Comm + Spatial MPI communicator. + + + Attributes + ---------- + domain : dolfinx.mesh.Mesh + Computational mesh. + W : FunctionSpace + Mixed function space for velocity and pressure. + V, Q : FunctionSpace + Velocity and pressure subspaces. + w : Function + Solution vector (velocity–pressure). + rhs : Function + Right-hand side of the nonlinear system. + F : UFL form + Nonlinear residual form. + M : PETSc.Mat + Mass matrix (velocity block). + Mf : PETSc.Mat + Mass matrix for the monolithic system (velocity and pressure). + + """ + dtype_u = mesh + dtype_f = mesh + + def __init__(self, nelems=128, t0=0.0, family='CG', order=4, nu=0.1, comm=MPI.COMM_WORLD): + + # Create mesh + domain = dfx.mesh.create_rectangle(comm, + [np.array([-0.5, -0.5]), np.array([0.5, 0.5])], + [nelems, nelems], + cell_type=dfx.mesh.CellType.triangle, + ) + + # Define mixed function space for velocity and pressure + Ve = basix.ufl.element(family, domain.topology.cell_name(), order, shape=(domain.geometry.dim,)) + Pe = basix.ufl.element(family, domain.topology.cell_name(), order - 1, shape=()) + We = basix.ufl.mixed_element([Ve, Pe]) + self.W = dfx.fem.functionspace(domain, We) + self.V, WV_map = self.W.sub(0).collapse() + self.Q, WQ_map = self.W.sub(1).collapse() + + # get number of dofs and set up the problem + self.tmp_u = dfx.fem.Function(self.W) + nx = len(self.tmp_u.x.array) + print('DoFs on this level:', nx) + + # invoke super init, passing number of dofs, dtype_u and dtype_f + super().__init__(init=(nx, comm, np.dtype('float64'))) + self._makeAttributeAndRegister('nelems', 't0', 'family', 'order', 'nu', 'comm', localVars=locals(), readOnly=True) + + # initialize functions for the solution, source term, and right-hand side + w_zero = dfx.fem.Function(self.W) + self.g = dfx.fem.Function(self.V) + self.rhs = dfx.fem.Function(self.W) + self.tmp_f = dfx.fem.Function(self.W) + self.u_bndry = dfx.fem.Function(self.W) + + # interpolate initial and boundary conditions + self.convert_to_fenicsx_vector(input=self.u_exact(t0), output=self.u_bndry) + + w_zero.sub(0).interpolate(lambda x: np.vstack(( np.full(x.shape[1], 0.0, dtype=np.float64),np.full(x.shape[1], 0.0, dtype=np.float64)))) + w_zero.sub(1).interpolate(lambda x: np.full(x.shape[1], 0.0, dtype=np.float64)) + + # locate boundary facets + fdim = domain.topology.dim - 1 + + w_dbc = dfx.fem.Function(self.W) + w_dbc.sub(0).interpolate(lambda x: np.vstack((np.full(x.shape[1], 1.0, dtype=np.float64),np.full(x.shape[1], 0.0, dtype=np.float64)))) + w_dbc.sub(1).interpolate(lambda x: np.full(x.shape[1], 1.0, dtype=np.float64)) + # + TopBottom = dfx.mesh.locate_entities_boundary(domain, fdim, lambda x: np.logical_or(np.isclose(x[1], -0.5) , np.isclose(x[1], 0.5))) + LeftRight = dfx.mesh.locate_entities_boundary(domain, fdim, lambda x: np.logical_or(np.isclose(x[0], -0.5) , np.isclose(x[0], 0.5))) + + dofs_TopBottom_u = dfx.fem.locate_dofs_topological((self.W.sub(0), self.V), fdim, TopBottom) + dofs_TopBottom_p = dfx.fem.locate_dofs_topological((self.W.sub(1), self.Q), fdim, TopBottom) + self.dofs_LeftRight_u = dfx.fem.locate_dofs_topological((self.W.sub(0), self.V), fdim, LeftRight) + self.dofs_LeftRight_p = dfx.fem.locate_dofs_topological((self.W.sub(1), self.Q), fdim, LeftRight) + # + bcu_TB = dfx.fem.dirichletbc(w_dbc.sub(0).collapse(), dofs_TopBottom_u, self.W.sub(0)) + bcp_TB = dfx.fem.dirichletbc(w_dbc.sub(1).collapse(), dofs_TopBottom_p, self.W.sub(1)) + self.bcs = [bcu_TB, bcp_TB] + + # homogeneous boundary conditions for the velocity and pressure for the residual + bndry_facets = dfx.mesh.locate_entities_boundary(domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool)) + dofs_bndry_u = dfx.fem.locate_dofs_topological((self.W.sub(0), self.V), fdim, bndry_facets) + dofs_bndry_p = dfx.fem.locate_dofs_topological((self.W.sub(1), self.Q), fdim, bndry_facets) + bcu_hom = dfx.fem.dirichletbc(w_zero.sub(0).collapse(), dofs_bndry_u, self.W.sub(0)) + bcp_hom = dfx.fem.dirichletbc(w_zero.sub(1).collapse(), dofs_bndry_p, self.W.sub(1)) + self.bc_hom = [bcu_hom, bcp_hom] + self.fix_bc_for_residual = True + + # Define trial and test functions + self.u, self.p = ufl.TrialFunctions(self.W) + self.v, self.q = ufl.TestFunctions(self.W) + + # Assemble mass matrices + ''' + | M_v 0 | | M_v 0 | + M = | | and Mf = | | + | 0 0 | | 0 M_p | + + where M_v and M_p are the mass matrix in the velocity and pressure spaces, respectively. + ''' + a_M = ufl.dot(self.u, self.v) * ufl.dx + self.M = dolfinx.fem.petsc.assemble_matrix(dfx.fem.form(a_M), bcs=[]) + self.M.assemble() + + a_Mf = ufl.dot(self.u, self.v) * ufl.dx + ufl.dot(self.p, self.q) * ufl.dx + self.Mf = dolfinx.fem.petsc.assemble_matrix(dfx.fem.form(a_Mf), bcs=[]) + self.Mf.assemble() + + self.Lin_solver = PETSc.KSP().create(domain.comm) + self.Lin_solver.setType(PETSc.KSP.Type.BCGS) + self.Lin_solver.setOperators(self.Mf) + self.Lin_solver.setTolerances(rtol=1e-10, atol=1e-10) + pc = self.Lin_solver.getPC() + pc.setType(PETSc.PC.Type.JACOBI) + + + # Initialize the nonlinear form for the Navier–Stokes equations + self.factor = dfx.fem.Constant(domain, PETSc.ScalarType(0.0)) + + self.w = dfx.fem.Function(self.W) + + u, p = ufl.split(self.w) + rhs_u, rhs_p = ufl.split(self.rhs) + + # define the nonlinear form to be compatible with the multi-point constraints + F = (1/self.factor) * ufl.dot(u, self.v) * ufl.dx + F += ufl.dot(ufl.dot(u, ufl.nabla_grad(u)), self.v) * ufl.dx + F += ufl.inner(self.nu * ufl.nabla_grad(u), ufl.nabla_grad(self.v)) * ufl.dx + F -= ufl.dot(p, ufl.div(self.v)) * ufl.dx + F -= ufl.dot(self.g, self.v) * ufl.dx + F -= ufl.dot(ufl.div(u), self.q) * ufl.dx + F -= (1/self.factor) * ufl.dot(rhs_u, self.v) * ufl.dx + F -= (1/self.factor) * ufl.dot(rhs_p, self.q) * ufl.dx + self.F = F + + # Assemble the nonlinear problem and configure the Newton solver + self.petsc_options = { + "snes_type": "newtonls", + "ksp_type": "preonly", + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps", + "snes_atol": 1e-12, + "snes_rtol": 1e-12, + "snes_linesearch_type": "none", + "ksp_error_if_not_converged": True, + "snes_error_if_not_converged": True, + } + + # plotter for solution visualization + self.plotter = None + + @staticmethod + def convert_to_fenicsx_vector(input, output): + output.x.array[:] = input[:] + + @staticmethod + def convert_from_fenicsx_vector(input, output): + output[:] = input.x.array[:] + + def solve_system(self, rhs, factor, u0, t): + r""" + Dolfin's nonlinear solver for : + .. math:: + (u,v) + factor * (u \cdot \nabla u, v) + factor * \nu (\nabla u, \nabla v) - factor * (p, \nabla \cdot v) - factor * (g, v) - factor * (div(u), q) = (rhs_u, v) + (rhs_p, q) + + Parameters + ---------- + rhs : dtype_f + Right-hand side for the nonlinear system. + factor : float + Abbrev. for the node-to-node stepsize (or any other factor required). + u0 : dtype_u + Initial guess for the iterative solver (not used here so far). + t : float + Current time. + + Returns + ------- + w : dtype_u + Solution. + """ + + w = self.dtype_u(self.init) + self.factor.value = factor + + # update the source term at time t + self.g.interpolate(self.source_term(t)) + + # get boundary values from the exact solution at time t and update the boundary conditions + self.convert_to_fenicsx_vector(input=self.u_exact(t), output=self.u_bndry) + # + bcu_LR = dfx.fem.dirichletbc(self.u_bndry.sub(0).collapse(), self.dofs_LeftRight_u, self.W.sub(0)) + bcp_LR = dfx.fem.dirichletbc(self.u_bndry.sub(1).collapse(), self.dofs_LeftRight_p, self.W.sub(1)) + bcup = self.bcs + [bcu_LR, bcp_LR] + # + self.convert_to_fenicsx_vector(input=rhs, output=self.tmp_f) + self.tmp_f = self.__invert_mass_matrix(self.tmp_f) + + self.rhs.x.array[:] = self.tmp_f.x.array[:] + + problem = dolfinx.fem.petsc.NonlinearProblem(self.F, self.w, bcs=bcup, + petsc_options=self.petsc_options, petsc_options_prefix="nonlinearsolver") + + # solve the nonlinear system using Dolfinx's nonlinear solver + problem.solve() + + self.w.x.scatter_forward() + + # + self.convert_from_fenicsx_vector(input=self.w, output=w) + + + return w + + def eval_f(self, w, t): + """ + Routine to evaluate both parts of the right-hand side of the problem. + + Parameters + ---------- + w : dtype_u + Current values of the numerical solution. + t : float + Current time at which the numerical solution is computed. + + Returns + ------- + f : dtype_f + The right-hand side. + """ + + f = self.dtype_f(self.init) + + self.convert_to_fenicsx_vector(input=w, output=self.tmp_u) + + u, p = ufl.split(self.tmp_u) + + # get the forcing term + g = self.source_term(t) + + F = -ufl.dot(ufl.dot(u, ufl.nabla_grad(u)), self.v) * ufl.dx + F -= ufl.inner(self.nu * ufl.nabla_grad(u), ufl.nabla_grad(self.v)) * ufl.dx + F += ufl.dot(p, ufl.div(self.v)) * ufl.dx + F += ufl.dot(g, self.v) * ufl.dx + F += ufl.dot(ufl.div(u), self.q) * ufl.dx + + b = dolfinx.fem.petsc.assemble_vector(dfx.fem.form(F)) + b.assemble() + b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + + # + with b.localForm() as loc_b, self.tmp_f.x.petsc_vec.localForm() as loc_f: + loc_b.copy(loc_f) + + self.convert_from_fenicsx_vector(input=self.tmp_f, output=f) + + return f + + def apply_mass_matrix(self, u): + r""" + Routine to apply mass matrix. + + Parameters + ---------- + u : dtype_u + Current values of the numerical solution. + + Returns + ------- + me : dtype_u + The product :math:`M \vec{u}`. + """ + + self.convert_to_fenicsx_vector(input=u, output=self.tmp_u) + self.M.mult(self.tmp_u.x.petsc_vec, self.tmp_f.x.petsc_vec) + me = self.dtype_u(self.init) + self.convert_from_fenicsx_vector(input=self.tmp_f, output=me) + return me + + def __invert_mass_matrix(self, u): + r""" + Helper routine to invert mass matrix. + + Parameters + ---------- + u : dtype_u + Current values of the numerical solution. + + Returns + ------- + me : dtype_u + The product :math:`Mf^{-1} \vec{u}`. + """ + self.Lin_solver.solve(u.x.petsc_vec, self.tmp_f.x.petsc_vec) + self.tmp_f.x.scatter_forward() + + return self.tmp_f + + + def source_term(self, t): + r""" + Routine to compute the source term at time :math:`t`. + + Parameters + ---------- + t : float + Time of the exact solution. + + Returns + ------- + g : dtype_u + Source term. + """ + + self.g.interpolate( + lambda x: np.vstack(( + np.pi / 34 * np.exp(-16 * np.pi**2 * self.nu * t) * np.sin(4 * np.pi * (t - x[0])) * \ + np.cos(np.pi * x[1]) * (32 - 17*np.cos(np.pi*x[1])), + # + 2 * np.pi**2 * self.nu * np.exp(-8 * np.pi**2 * self.nu * t) * np.cos(2 * np.pi * (t - x[0])) \ + - np.pi * np.exp(-16 * np.pi**2 * self.nu * t) * np.sin(np.pi * x[1]) * (2 * np.cos(np.pi \ + * x[1])**3 + 4/17 * np.cos(4 * np.pi * (t - x[0]))) + )) + ) + return self.g + + + def u_exact(self, t): + r""" + Routine to compute the exact solution at time :math:`t`. + + Parameters + ---------- + t : float + Time of the exact solution. + + Returns + ------- + me : dtype_u + Exact solution. + """ + + + w = dfx.fem.Function(self.W) + + w.sub(0).interpolate(lambda x: np.vstack(( + 1.0 - np.exp(-8 * np.pi**2 * self.nu * t) * np.sin(2.0 * np.pi * (x[0] - t)) \ + * np.sin(np.pi * x[1]) * np.cos(np.pi * x[1]), + # + - np.exp(-8 * np.pi**2 * self.nu * t) * np.cos(2.0 * np.pi * (x[0] - t)) \ + * (np.cos(np.pi * x[1]))**2 + )) + ) + + w.sub(1).interpolate( lambda x: (4/17) * np.exp(-16 * np.pi**2 * self.nu * t) * \ + np.cos(4 * np.pi * (x[0]-t)) * np.cos(np.pi * x[1]) + 1 + ) + + + + me = self.dtype_u(self.init) + self.convert_from_fenicsx_vector(input=w, output=me) + + return me + + def fix_residual(self, res): + """ + Applies homogeneous Dirichlet boundary conditions to the residual + + Parameters + ---------- + res : dtype_u + Residual + """ + self.convert_to_fenicsx_vector(input=res, output=self.tmp_u) + dolfinx.fem.petsc.set_bc(self.tmp_u.x.petsc_vec, self.bc_hom) + self.tmp_u.x.scatter_forward() + self.convert_from_fenicsx_vector(input=self.tmp_u, output=res) + return None + + +class fenicsx_NSE_periodic_mass(fenicsx_NSE_mass): + r""" + Monolithic incompressible Navier–Stokes problem in 2D using FEniCSx. + + This class implements a finite element discretization of the incompressible + Navier–Stokes equations in a square domain using a monolithic formulation + for velocity and pressure. + + Governing equations: + .. math:: + \frac{\partial \mathbf{u}}{\partial t} + + (\mathbf{u} \cdot \nabla)\mathbf{u} + - \nu \Delta \mathbf{u} + + \nabla p = f, + \quad \nabla \cdot \mathbf{u} = 0 + + Domain: + A square domain :math:`\Omega = [-0.5, 0.5]^2` discretized with triangular + elements. + + Spatial discretization: + - Continuous Lagrange elements for velocity and pressure + - Mixed finite element space :math:`W = V \times Q` + - Monolithic formulation for the coupled velocity–pressure system + + Manufactured solution: + A smooth, analytical, time-dependent solution is given by: + ..math:: + u(x,y,t) = 1 - e^{-8\pi^{2}\nu t}\sin(2\pi(x-t))\sin(\pi y)\cos(\pi y) + v(x,y,t) = - e^{-8\pi^{2}\nu t} \cos\(2\pi(x-t)) \cos^{2}(\pi y) + p(x,y,t) = 1 + \frac{4}{17} e^{-16\pi^{2}\nu t} \cos\(4\pi(x-t)) \cos(\pi y) + + The solution is specifically constructed to investigate the impact of boundary + conditions on numerical accuracy. The source term g =(g_1, g_2) is derived from + the exact solution to ensure that it satisfies the Navier–Stokes equations. + + Boundary conditions: + - periodic boundary conditions on the left and right boundaries + - Time-independent Dirichlet conditions on the top and bottom boundaries + + In this class, the same problem as "fenicsx_NSE_mass" is formulated with periodic + boundary conditions on the left and right boundaries, while retaining the time- + independent Dirichlet conditions on the top and bottom boundaries, in order to + isolate and compare the effects of time-dependent boundary data on the numerical + accuracy. + + + Parameters + ---------- + nelems : int + Number of elements per spatial direction. + t0 : float + Initial time. + family : str + Finite element family (e.g., 'CG' for Lagrange). + order : int + Polynomial degree of the function spaces. + nu : float + Kinematic viscosity. + comm : MPI.Comm + Spatial MPI communicator. + + + Attributes + ---------- + domain : dolfinx.mesh.Mesh + Computational mesh. + W : FunctionSpace + Mixed function space for velocity and pressure. + V, Q : FunctionSpace + Velocity and pressure subspaces. + w : Function + Solution vector (velocity–pressure). + rhs : Function + Right-hand side of the nonlinear system. + F : UFL form + Nonlinear residual form. + M : PETSc.Mat + Mass matrix (velocity block). + Mf : PETSc.Mat + Mass matrix for the monolithic system (velocity and pressure). + """ + + def __init__(self, nelems=128, t0=0.0, family='CG', order=4, nu=0.1, comm=MPI.COMM_WORLD): + super().__init__(nelems=nelems, t0=t0, family=family, order=order, nu=nu, comm=comm) + + # get the mesh and its geometric properties for periodicity + domain = self.W.mesh + + self.xmin = domain.geometry.x[:, 0].min() + self.xmax = domain.geometry.x[:, 0].max() + self.ymin = domain.geometry.x[:, 1].min() + self.ymax = domain.geometry.x[:, 1].max() + self.Lx = float(self.xmax - self.xmin) + self.Ly = float(self.ymax - self.ymin) + + fdim = domain.topology.dim - 1 + + # locate facets on the right boundary for periodicity and create meshtags for them + facets_x = dfx.mesh.locate_entities_boundary(domain, fdim, self.periodic_boundary_x) + arg_sort_x = np.argsort(facets_x) + mt_x = dfx.mesh.meshtags(domain, fdim, facets_x[arg_sort_x], np.full(len(facets_x), 2, dtype=np.int32)) + + # create multi-point constraints for periodicity and add the periodic constraints to tthe boundary conditions + self.mpc = dolfinx_mpc.MultiPointConstraint(self.W) + self.mpc.create_periodic_constraint_topological(self.W.sub(0).sub(0), mt_x, 2, self.periodic_relation_x, self.bcs) + self.mpc.create_periodic_constraint_topological(self.W.sub(0).sub(1), mt_x, 2, self.periodic_relation_x, self.bcs) + self.mpc.create_periodic_constraint_topological(self.W.sub(1), mt_x, 2, self.periodic_relation_x, self.bcs) + self.mpc.finalize() + + # reassemble the nonlinear problem + self.factor = dfx.fem.Constant(domain, PETSc.ScalarType(0.0)) + + # redefine the solution function to be compatible with the multi-point constraints + W_mpc = self.mpc.function_space + self.w = dfx.fem.Function(W_mpc) + u, p = ufl.split(self.w) + rhs_u, rhs_p = ufl.split(self.rhs) + + # redefine the nonlinear form to be compatible with the multi-point constraints + F = (1/self.factor) * ufl.dot(u, self.v) * ufl.dx + F += ufl.dot(ufl.dot(u, ufl.nabla_grad(u)), self.v) * ufl.dx + F += ufl.inner(self.nu * ufl.nabla_grad(u), ufl.nabla_grad(self.v)) * ufl.dx + F -= ufl.dot(p, ufl.div(self.v)) * ufl.dx + F -= ufl.dot(self.g, self.v) * ufl.dx + F -= ufl.dot(ufl.div(u), self.q) * ufl.dx + F -= (1/self.factor) * ufl.dot(rhs_u, self.v) * ufl.dx + F -= (1/self.factor) * ufl.dot(rhs_p, self.q) * ufl.dx + + # compute the Jacobian of the nonlinear form for use in the Newton solver + Jac = ufl.derivative(F, self.w, ufl.TrialFunction(self.W)) + + # assemble the nonlinear problem with the multi-point constraints and configure the Newton solver + self.problem = dolfinx_mpc.NonlinearProblem(F, self.w, mpc=self.mpc, bcs=self.bcs, J=Jac, petsc_options=self.petsc_options) + + def solve_system(self, rhs, factor, u0, t): + r""" + Dolfin's nonlinear solver for : + .. math:: + (u,v) + factor * (u \cdot \nabla u, v) + factor * \nu (\nabla u, \nabla v) - factor * (p, \nabla \cdot v) - factor * (g, v) - factor * (div(u), q) = (rhs_u, v) + (rhs_p, q) + + Parameters + ---------- + rhs : dtype_f + Right-hand side for the nonlinear system. + factor : float + Abbrev. for the node-to-node stepsize (or any other factor required). + u0 : dtype_u + Initial guess for the iterative solver (not used here so far). + t : float + Current time. + + Returns + ------- + w : dtype_u + Solution. + """ + self.g.interpolate(self.source_term(t)) + + u = self.dtype_u(self.init) + self.factor.value = factor + # + self.convert_to_fenicsx_vector(input=rhs, output=self.tmp_f) + self.tmp_f = super()._fenicsx_NSE_mass__invert_mass_matrix(self.tmp_f) + # + self.rhs.x.array[:] = self.tmp_f.x.array[:] + # + self.problem.solve() + self.w.x.scatter_forward() + # + self.convert_from_fenicsx_vector(input=self.w, output=u) + + return u + + def periodic_boundary_x(self, x): + r""" + Identifies the left and right boundaries for periodicity. + Parameters + ---------- + x : array_like + Coordinates of the points to check for periodicity. + Returns + ------- + boolTrue if the point is on the right boundary (x = xmax), False otherwise. + + """ + eps = 10 * np.finfo(x.dtype).resolution + return np.isclose(x[0], self.xmax, atol=eps) + + def periodic_relation_x(self, x): + r""" + Maps points from the right boundary to the left boundary for periodicity. + Parameters + ---------- + x : array_like + Coordinates of the points on the right boundary to be mapped. + Returns + ------- + out_x : array_like + Coordinates of the corresponding points on the left boundary after applying periodicity. + """ + out_x = np.zeros_like(x) + out_x[0] = x[0] - self.Lx + out_x[1] = x[1] + out_x[2] = x[2] + return out_x + + \ No newline at end of file diff --git a/pySDC/projects/StroemungsRaum_FEniCSx/run_Navier_Stokes_equations_TaylGreen_monolithic_FEniCSx.py b/pySDC/projects/StroemungsRaum_FEniCSx/run_Navier_Stokes_equations_TaylGreen_monolithic_FEniCSx.py new file mode 100755 index 0000000000..27f879f7a2 --- /dev/null +++ b/pySDC/projects/StroemungsRaum_FEniCSx/run_Navier_Stokes_equations_TaylGreen_monolithic_FEniCSx.py @@ -0,0 +1,169 @@ +import sys +import numpy as np +import dolfinx as dfx +import ufl +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI +from pySDC.projects.StroemungsRaum.sweepers.generic_implicit_mass import generic_implicit_mass +from pySDC.projects.StroemungsRaum_FEniCSx.problem_classes.NavierStokes_2D_TaylGreen_monolithic_FEniCSx import ( + fenicsx_NSE_mass, + fenicsx_NSE_periodic_mass + ) + +def setup(t0=0, periodic=False): + """ + Helper routine to set up parameters + + Args: + t0: float, + initial time + periodic: bool, + whether to use periodic boundary conditions or not + + Returns: + description: dict, + pySDC description dictionary containing problem and method parameters. + controller_params: dict, + Parameters for the pySDC controller. + """ + # time step size + dt = 0.1 + + # initialize level parameters + level_params = dict() + level_params['restol'] = 1e-10 + level_params['dt'] = dt + + # initialize step parameters + step_params = dict() + step_params['maxiter'] = 15 + + # initialize sweeper parameters + sweeper_params = dict() + sweeper_params['quad_type'] = 'RADAU-RIGHT' + sweeper_params['num_nodes'] = [2] + sweeper_params['QI'] = ['LU'] + + problem_params = dict() + problem_params['nu'] = 0.2 + problem_params['t0'] = t0 + problem_params['nelems'] = [64] + problem_params['family'] = 'CG' + problem_params['order'] = [2] + + # initialize controller parameters + controller_params = dict() + controller_params['logger_level'] = 20 + + # Fill description dictionary for easy hierarchy creation + description = dict() + if periodic: + description['problem_class'] = fenicsx_NSE_periodic_mass + else: + description['problem_class'] = fenicsx_NSE_mass + description['sweeper_class'] = generic_implicit_mass + description['problem_params'] = problem_params + description['sweeper_params'] = sweeper_params + description['level_params'] = level_params + description['step_params'] = step_params + + return description, controller_params + + +def run_simulation(description, controller_params, Tend): + """ + Run the time integration for the 2D Navier-Stokes equations. + + Args: + description: dict, + pySDC problem and method description. + controller_params: dict, + Parameters for the pySDC controller. + Tend: float, + Final simulation time. + + Returns: + P: problem instance, + Problem instance containing the final solution and other problem-related information. + stats: dict, + collected runtime statistics, + uend: FEniCS function, + Final solution at time Tend. + """ + # get initial time from description + t0 = description['problem_params']['t0'] + + # quickly generate block of steps + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + + # get initial values on finest level + P = controller.MS[0].levels[0].prob + uinit = P.u_exact(t0) + + # call main function to get things done... + uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + + return P, stats, uend + +def run_postprocessing(P, uend, Tend): + """ + Postprocess and store simulation results for visualization and analysis. + + Args: + P: Problem instance, + Problem instance containing the final solution and other problem-related information. + uend: FEniCS function, + Final solution at time Tend. + Tend: float, + Final simulation time. + + Returns: + rel_error_u: float, + Relative L2 error in velocity compared to the exact solution. + rel_error_p: float, + Relative L2 error in pressure compared to the exact solution. + """ + wx = dfx.fem.Function(P.W) + wn = dfx.fem.Function(P.W) + # + wx.x.array[:] = P.u_exact(Tend)[:] + wn.x.array[:] = uend[:] + # + un, pn = wn.split() + un_ = un.collapse() + pn_ = pn.collapse() + # + ue, pe = wx.split() + ue_ = ue.collapse() + pe_ = pe.collapse() + # + error_u = np.sqrt(dfx.fem.assemble_scalar(dfx.fem.form(ufl.dot(un_ - ue_, un_ - ue_) * ufl.dx))) + error_p = np.sqrt(dfx.fem.assemble_scalar(dfx.fem.form((pn_ - pe_) * (pn_ - pe_) * ufl.dx))) + + # + norm_u = np.sqrt(dfx.fem.assemble_scalar(dfx.fem.form(ufl.dot(ue_, ue_) * ufl.dx))) + norm_p = np.sqrt(dfx.fem.assemble_scalar(dfx.fem.form(pe_ * pe_ * ufl.dx))) + + # + rel_error_u = error_u / norm_u + rel_error_p = error_p / norm_p + + # + print(f"L2 error in velocity : {error_u:.3e} (relative: {rel_error_u:.3e})") + print(f"L2 error in pressure : {error_p:.3e} (relative: {rel_error_p:.3e})") + + return rel_error_u, rel_error_p + + +if __name__ == "__main__": + + t0 = 0.0 + Tend =0.1 + + # run the setup to get description and controller parameters + description, controller_params = setup(t0=t0, periodic=True) + + # run the simulation and get the problem, stats and final solution + problem, stats, uend = run_simulation(description, controller_params, Tend) + + # run postprocessing to save parameters and solutions for visualization + rel_error_u, rel_error_p = run_postprocessing(problem, uend, Tend) \ No newline at end of file diff --git a/pySDC/projects/StroemungsRaum_FEniCSx/tests/test_Navier_Stokes_TaylGreen_monolithic_FEniCSx.py b/pySDC/projects/StroemungsRaum_FEniCSx/tests/test_Navier_Stokes_TaylGreen_monolithic_FEniCSx.py new file mode 100644 index 0000000000..a40ef04525 --- /dev/null +++ b/pySDC/projects/StroemungsRaum_FEniCSx/tests/test_Navier_Stokes_TaylGreen_monolithic_FEniCSx.py @@ -0,0 +1,154 @@ +import pytest +import numpy as np +import ufl +import dolfinx as dfx +import dolfinx.fem.petsc + +from pySDC.projects.StroemungsRaum_FEniCSx.problem_classes.NavierStokes_2D_TaylGreen_monolithic_FEniCSx import ( + fenicsx_NSE_mass, + fenicsx_NSE_periodic_mass, + ) + +@pytest.mark.fenicsx +def test_solve_system(): +# + # Physical and numerical parameters for the test case + nu = 0.02 + t = 0.0 + factor = 0.01 + + # Create the 2D monolithic Navier-Stokes problem for both non-periodic and periodic cases + prob = fenicsx_NSE_mass(t0=0.0, nelems=32, family='CG', order=2, nu=nu) + prob_periodic = fenicsx_NSE_periodic_mass(t0=0.0, nelems=32, family='CG', order=2, nu=nu) + + # Evaluate the source term at the given time + g = prob.source_term(t) + + # get the exact solution at time t + wex = dfx.fem.Function(prob.W) + wex.x.array[:] = prob.u_exact(t)[:] + + # Split the exact solution into velocity and pressure components + u, p = ufl.split(wex) + + # assemble the right-hand side vector + rhs_f = ufl.dot(u, prob.v) * ufl.dx + rhs_f += factor * ufl.dot(ufl.dot(u, ufl.nabla_grad(u)), prob.v) * ufl.dx + rhs_f += factor * ufl.inner(prob.nu * ufl.nabla_grad(u), ufl.nabla_grad(prob.v)) * ufl.dx + rhs_f -= factor * ufl.dot(p, ufl.div(prob.v)) * ufl.dx + rhs_f -= factor * ufl.dot(g, prob.v) * ufl.dx + rhs_f -= factor * ufl.dot(ufl.div(u), prob.q) * ufl.dx + + b = dolfinx.fem.petsc.assemble_vector(dfx.fem.form(rhs_f)) + b.assemble() + #b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + + + # Solve the linear system for both non-periodic and periodic cases + w = prob.solve_system(rhs=b, factor=factor, u0=prob.dtype_u(prob.init), t=t) + w_periodic = prob_periodic.solve_system(rhs=b, factor=factor, u0=prob.dtype_u(prob.init), t=t) + + # get the exact solution at time t for error computation + wx = prob.dtype_u(prob.init) + wx[:] = wex.x.array[:] + + # Compute the relative monolithic velocity–pressure error for both non-periodic and periodic cases + rel_err = abs(wx - w) / abs(wx) + rel_err_periodic = abs(wx - w_periodic) / abs(wx) + assert rel_err < 1e-5, (f"Relative monolithic velocity–pressure error {rel_err:.3e} exceeds tolerance") + assert rel_err_periodic < 1e-5, ( + f"Relative monolithic velocity–pressure error (periodic) {rel_err_periodic:.3e} exceeds tolerance") + + +@pytest.mark.fenicsx +def test_eval_f(): + + # parameters for the test case + nu = 0.2 + t = 0.0 + + # create the 2D monolithic Navier-Stokes problem + prob = fenicsx_NSE_mass(t0=0.0, nelems=64, family='CG', order=2, nu=nu) + + # get the exact solution at time t + wex = dfx.fem.Function(prob.W) + wex.sub(0).interpolate(lambda x: + np.vstack((-1.0 *np.cos(x[0])*np.sin(x[1])*np.exp(-2.0*nu*t),np.sin(x[0])*np.cos(x[1])*np.exp(-2.0*nu*t)))) + wex.sub(1).interpolate( lambda x: -0.25*(np.cos(2.0*x[0]) + np.cos(2.0*x[1]))*np.exp(-4.0*nu*t)) + + # evaluate the right-hand side vector using the problem's eval_f method + f = prob.eval_f(w=wex.x.array[:], t=t) + + # compute the expected right-hand side vector by applying the mass matrix to the exact solution + fw = dfx.fem.Function(prob.W) + fw.sub(0).interpolate(lambda x: np.vstack(( + 2.0*nu*np.cos(x[0])*np.sin(x[1])*np.exp(-2.0*nu*t) + + np.pi / 34 * np.exp(-16 * np.pi**2 * nu * t) * np.sin(4 * np.pi * (t - x[0])) + * np.cos(np.pi * x[1]) * (32 - 17*np.cos(np.pi*x[1])), + # + -2.0*nu*np.sin(x[0])*np.cos(x[1])*np.exp(-2.0*nu*t) + + 2 * np.pi**2 * nu * np.exp(-8 * np.pi**2 * nu * t) * np.cos(2 * np.pi * (t - x[0])) + - np.pi * np.exp(-16 * np.pi**2 * nu * t) * np.sin(np.pi * x[1]) * (2 * np.cos(np.pi \ + * x[1])**3 + 4/17 * np.cos(4 * np.pi * (t - x[0]))) + ))) + fw.sub(1).interpolate(lambda x: np.full(x.shape[1], 0.0, dtype=np.float64)) + + # apply the mass matrix to obtain the expected right-hand side vector + Mfw = dfx.fem.Function(prob.W) + prob.M.mult(fw.x.petsc_vec, Mfw.x.petsc_vec) + + # convert the expected right-hand side vector to a numpy array for comparison + fwx = prob.dtype_u(prob.init) + fwx[:] = Mfw.x.array[:] + + # apply boundary conditions to both vectors before comparison + prob.fix_residual(f) + prob.fix_residual(fwx) + + # compute the relative error between computed and expected right-hand sides + rel_err = abs(f - fwx) / abs(fwx) + print(f"Relative error in right-hand side evaluation: {rel_err:.3e}") + assert rel_err < 1e-5, f"Relative error {rel_err} exceeds tolerance" + + +@pytest.mark.fenicsx +def test_problem_class(): + + from pySDC.projects.StroemungsRaum_FEniCSx.run_Navier_Stokes_equations_TaylGreen_monolithic_FEniCSx import ( + setup, + run_simulation, + run_postprocessing + ) + + t0 = 0.0 + Tend = 0.1 + + # -- Non-periodic case -- + # run the setup to get description and controller parameters + description, controller_params = setup(t0=t0, periodic=False) + + # run the simulation and get the problem, stats and final solution + problem, stats, uend = run_simulation(description, controller_params, Tend) + + # run postprocessing to save parameters and solutions for visualization + rel_error_u, rel_error_p = run_postprocessing(problem, uend, Tend) + + print(f"Relative L2 error in velocity: {rel_error_u:.3e}") + print(f"Relative L2 error in pressure: {rel_error_p:.3e}") + # assert that the relative errors in velocity and pressure are within acceptable limits + assert rel_error_u < 1e-2, f"Error in velocity error {rel_error_u} exceeds tolerance" + assert rel_error_p < 6e-2, f"Error in pressure error {rel_error_p} exceeds tolerance" + + # -- Periodic case -- + # run the setup to get description and controller parameters + description, controller_params = setup(t0=t0, periodic=True) + + # run the simulation and get the problem, stats and final solution + problem, stats, uend = run_simulation(description, controller_params, Tend) + + # run postprocessing to save parameters and solutions for visualization + rel_error_u, rel_error_p = run_postprocessing(problem, uend, Tend) + + # assert that the relative errors in velocity and pressure are within acceptable limits + assert rel_error_u < 1e-2, f"Error in velocity error (periodic bc) {rel_error_u} exceeds tolerance" + assert rel_error_p < 1e-2, f"Error in pressure error (periodic bc) {rel_error_p} exceeds tolerance" \ No newline at end of file From f2163b42a8ae17af59bb738c812bf0b618216d4f Mon Sep 17 00:00:00 2001 From: Ouardghi Date: Mon, 18 May 2026 17:13:28 +0200 Subject: [PATCH 2/2] fix linting --- ...rStokes_2D_TaylGreen_monolithic_FEniCSx.py | 229 ++++++++++-------- ..._equations_TaylGreen_monolithic_FEniCSx.py | 42 ++-- ...ier_Stokes_TaylGreen_monolithic_FEniCSx.py | 91 ++++--- 3 files changed, 207 insertions(+), 155 deletions(-) diff --git a/pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py b/pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py index 9149e4bcb9..8b249fd0fb 100755 --- a/pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py +++ b/pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py @@ -12,7 +12,7 @@ class fenicsx_NSE_mass(Problem): - r""" + r""" Monolithic incompressible Navier–Stokes problem in 2D using FEniCSx. This class implements a finite element discretization of the incompressible @@ -41,10 +41,10 @@ class fenicsx_NSE_mass(Problem): u(x,y,t) = 1 - e^{-8\pi^{2}\nu t}\sin(2\pi(x-t))\sin(\pi y)\cos(\pi y) v(x,y,t) = - e^{-8\pi^{2}\nu t} \cos\(2\pi(x-t)) \cos^{2}(\pi y) p(x,y,t) = 1 + \frac{4}{17} e^{-16\pi^{2}\nu t} \cos\(4\pi(x-t)) \cos(\pi y) - - The solution is specifically constructed to investigate the impact of boundary - conditions on numerical accuracy. The source term g =(g_1, g_2) is derived from - the exact solution to ensure that it satisfies the Navier–Stokes equations. + + The solution is specifically constructed to investigate the impact of boundary + conditions on numerical accuracy. The source term g =(g_1, g_2) is derived from + the exact solution to ensure that it satisfies the Navier–Stokes equations. Boundary conditions: - Time-dependent Dirichlet conditions on the left and right boundaries @@ -95,18 +95,20 @@ class fenicsx_NSE_mass(Problem): Mass matrix for the monolithic system (velocity and pressure). """ + dtype_u = mesh dtype_f = mesh - + def __init__(self, nelems=128, t0=0.0, family='CG', order=4, nu=0.1, comm=MPI.COMM_WORLD): - - # Create mesh - domain = dfx.mesh.create_rectangle(comm, - [np.array([-0.5, -0.5]), np.array([0.5, 0.5])], - [nelems, nelems], - cell_type=dfx.mesh.CellType.triangle, - ) - + + # Create mesh + domain = dfx.mesh.create_rectangle( + comm, + [np.array([-0.5, -0.5]), np.array([0.5, 0.5])], + [nelems, nelems], + cell_type=dfx.mesh.CellType.triangle, + ) + # Define mixed function space for velocity and pressure Ve = basix.ufl.element(family, domain.topology.cell_name(), order, shape=(domain.geometry.dim,)) Pe = basix.ufl.element(family, domain.topology.cell_name(), order - 1, shape=()) @@ -122,7 +124,9 @@ def __init__(self, nelems=128, t0=0.0, family='CG', order=4, nu=0.1, comm=MPI.CO # invoke super init, passing number of dofs, dtype_u and dtype_f super().__init__(init=(nx, comm, np.dtype('float64'))) - self._makeAttributeAndRegister('nelems', 't0', 'family', 'order', 'nu', 'comm', localVars=locals(), readOnly=True) + self._makeAttributeAndRegister( + 'nelems', 't0', 'family', 'order', 'nu', 'comm', localVars=locals(), readOnly=True + ) # initialize functions for the solution, source term, and right-hand side w_zero = dfx.fem.Function(self.W) @@ -130,22 +134,34 @@ def __init__(self, nelems=128, t0=0.0, family='CG', order=4, nu=0.1, comm=MPI.CO self.rhs = dfx.fem.Function(self.W) self.tmp_f = dfx.fem.Function(self.W) self.u_bndry = dfx.fem.Function(self.W) - + # interpolate initial and boundary conditions self.convert_to_fenicsx_vector(input=self.u_exact(t0), output=self.u_bndry) - w_zero.sub(0).interpolate(lambda x: np.vstack(( np.full(x.shape[1], 0.0, dtype=np.float64),np.full(x.shape[1], 0.0, dtype=np.float64)))) + w_zero.sub(0).interpolate( + lambda x: np.vstack( + (np.full(x.shape[1], 0.0, dtype=np.float64), np.full(x.shape[1], 0.0, dtype=np.float64)) + ) + ) w_zero.sub(1).interpolate(lambda x: np.full(x.shape[1], 0.0, dtype=np.float64)) - + # locate boundary facets - fdim = domain.topology.dim - 1 - + fdim = domain.topology.dim - 1 + w_dbc = dfx.fem.Function(self.W) - w_dbc.sub(0).interpolate(lambda x: np.vstack((np.full(x.shape[1], 1.0, dtype=np.float64),np.full(x.shape[1], 0.0, dtype=np.float64)))) + w_dbc.sub(0).interpolate( + lambda x: np.vstack( + (np.full(x.shape[1], 1.0, dtype=np.float64), np.full(x.shape[1], 0.0, dtype=np.float64)) + ) + ) w_dbc.sub(1).interpolate(lambda x: np.full(x.shape[1], 1.0, dtype=np.float64)) # - TopBottom = dfx.mesh.locate_entities_boundary(domain, fdim, lambda x: np.logical_or(np.isclose(x[1], -0.5) , np.isclose(x[1], 0.5))) - LeftRight = dfx.mesh.locate_entities_boundary(domain, fdim, lambda x: np.logical_or(np.isclose(x[0], -0.5) , np.isclose(x[0], 0.5))) + TopBottom = dfx.mesh.locate_entities_boundary( + domain, fdim, lambda x: np.logical_or(np.isclose(x[1], -0.5), np.isclose(x[1], 0.5)) + ) + LeftRight = dfx.mesh.locate_entities_boundary( + domain, fdim, lambda x: np.logical_or(np.isclose(x[0], -0.5), np.isclose(x[0], 0.5)) + ) dofs_TopBottom_u = dfx.fem.locate_dofs_topological((self.W.sub(0), self.V), fdim, TopBottom) dofs_TopBottom_p = dfx.fem.locate_dofs_topological((self.W.sub(1), self.Q), fdim, TopBottom) @@ -164,7 +180,7 @@ def __init__(self, nelems=128, t0=0.0, family='CG', order=4, nu=0.1, comm=MPI.CO bcp_hom = dfx.fem.dirichletbc(w_zero.sub(1).collapse(), dofs_bndry_p, self.W.sub(1)) self.bc_hom = [bcu_hom, bcp_hom] self.fix_bc_for_residual = True - + # Define trial and test functions self.u, self.p = ufl.TrialFunctions(self.W) self.v, self.q = ufl.TestFunctions(self.W) @@ -184,7 +200,7 @@ def __init__(self, nelems=128, t0=0.0, family='CG', order=4, nu=0.1, comm=MPI.CO a_Mf = ufl.dot(self.u, self.v) * ufl.dx + ufl.dot(self.p, self.q) * ufl.dx self.Mf = dolfinx.fem.petsc.assemble_matrix(dfx.fem.form(a_Mf), bcs=[]) self.Mf.assemble() - + self.Lin_solver = PETSc.KSP().create(domain.comm) self.Lin_solver.setType(PETSc.KSP.Type.BCGS) self.Lin_solver.setOperators(self.Mf) @@ -192,24 +208,23 @@ def __init__(self, nelems=128, t0=0.0, family='CG', order=4, nu=0.1, comm=MPI.CO pc = self.Lin_solver.getPC() pc.setType(PETSc.PC.Type.JACOBI) - # Initialize the nonlinear form for the Navier–Stokes equations self.factor = dfx.fem.Constant(domain, PETSc.ScalarType(0.0)) self.w = dfx.fem.Function(self.W) - + u, p = ufl.split(self.w) rhs_u, rhs_p = ufl.split(self.rhs) - + # define the nonlinear form to be compatible with the multi-point constraints - F = (1/self.factor) * ufl.dot(u, self.v) * ufl.dx + F = (1 / self.factor) * ufl.dot(u, self.v) * ufl.dx F += ufl.dot(ufl.dot(u, ufl.nabla_grad(u)), self.v) * ufl.dx F += ufl.inner(self.nu * ufl.nabla_grad(u), ufl.nabla_grad(self.v)) * ufl.dx F -= ufl.dot(p, ufl.div(self.v)) * ufl.dx F -= ufl.dot(self.g, self.v) * ufl.dx F -= ufl.dot(ufl.div(u), self.q) * ufl.dx - F -= (1/self.factor) * ufl.dot(rhs_u, self.v) * ufl.dx - F -= (1/self.factor) * ufl.dot(rhs_p, self.q) * ufl.dx + F -= (1 / self.factor) * ufl.dot(rhs_u, self.v) * ufl.dx + F -= (1 / self.factor) * ufl.dot(rhs_p, self.q) * ufl.dx self.F = F # Assemble the nonlinear problem and configure the Newton solver @@ -224,7 +239,7 @@ def __init__(self, nelems=128, t0=0.0, family='CG', order=4, nu=0.1, comm=MPI.CO "ksp_error_if_not_converged": True, "snes_error_if_not_converged": True, } - + # plotter for solution visualization self.plotter = None @@ -235,7 +250,7 @@ def convert_to_fenicsx_vector(input, output): @staticmethod def convert_from_fenicsx_vector(input, output): output[:] = input.x.array[:] - + def solve_system(self, rhs, factor, u0, t): r""" Dolfin's nonlinear solver for : @@ -276,18 +291,18 @@ def solve_system(self, rhs, factor, u0, t): self.tmp_f = self.__invert_mass_matrix(self.tmp_f) self.rhs.x.array[:] = self.tmp_f.x.array[:] - - problem = dolfinx.fem.petsc.NonlinearProblem(self.F, self.w, bcs=bcup, - petsc_options=self.petsc_options, petsc_options_prefix="nonlinearsolver") - + + problem = dolfinx.fem.petsc.NonlinearProblem( + self.F, self.w, bcs=bcup, petsc_options=self.petsc_options, petsc_options_prefix="nonlinearsolver" + ) + # solve the nonlinear system using Dolfinx's nonlinear solver problem.solve() - + self.w.x.scatter_forward() # self.convert_from_fenicsx_vector(input=self.w, output=w) - return w @@ -326,7 +341,7 @@ def eval_f(self, w, t): b = dolfinx.fem.petsc.assemble_vector(dfx.fem.form(F)) b.assemble() b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) - + # with b.localForm() as loc_b, self.tmp_f.x.petsc_vec.localForm() as loc_f: loc_b.copy(loc_f) @@ -334,7 +349,7 @@ def eval_f(self, w, t): self.convert_from_fenicsx_vector(input=self.tmp_f, output=f) return f - + def apply_mass_matrix(self, u): r""" Routine to apply mass matrix. @@ -355,7 +370,7 @@ def apply_mass_matrix(self, u): me = self.dtype_u(self.init) self.convert_from_fenicsx_vector(input=self.tmp_f, output=me) return me - + def __invert_mass_matrix(self, u): r""" Helper routine to invert mass matrix. @@ -374,8 +389,7 @@ def __invert_mass_matrix(self, u): self.tmp_f.x.scatter_forward() return self.tmp_f - - + def source_term(self, t): r""" Routine to compute the source term at time :math:`t`. @@ -392,18 +406,25 @@ def source_term(self, t): """ self.g.interpolate( - lambda x: np.vstack(( - np.pi / 34 * np.exp(-16 * np.pi**2 * self.nu * t) * np.sin(4 * np.pi * (t - x[0])) * \ - np.cos(np.pi * x[1]) * (32 - 17*np.cos(np.pi*x[1])), + lambda x: np.vstack( + ( + np.pi + / 34 + * np.exp(-16 * np.pi**2 * self.nu * t) + * np.sin(4 * np.pi * (t - x[0])) + * np.cos(np.pi * x[1]) + * (32 - 17 * np.cos(np.pi * x[1])), # - 2 * np.pi**2 * self.nu * np.exp(-8 * np.pi**2 * self.nu * t) * np.cos(2 * np.pi * (t - x[0])) \ - - np.pi * np.exp(-16 * np.pi**2 * self.nu * t) * np.sin(np.pi * x[1]) * (2 * np.cos(np.pi \ - * x[1])**3 + 4/17 * np.cos(4 * np.pi * (t - x[0]))) - )) + 2 * np.pi**2 * self.nu * np.exp(-8 * np.pi**2 * self.nu * t) * np.cos(2 * np.pi * (t - x[0])) + - np.pi + * np.exp(-16 * np.pi**2 * self.nu * t) + * np.sin(np.pi * x[1]) + * (2 * np.cos(np.pi * x[1]) ** 3 + 4 / 17 * np.cos(4 * np.pi * (t - x[0]))), + ) + ) ) return self.g - def u_exact(self, t): r""" Routine to compute the exact solution at time :math:`t`. @@ -418,24 +439,32 @@ def u_exact(self, t): me : dtype_u Exact solution. """ - - + w = dfx.fem.Function(self.W) - w.sub(0).interpolate(lambda x: np.vstack(( - 1.0 - np.exp(-8 * np.pi**2 * self.nu * t) * np.sin(2.0 * np.pi * (x[0] - t)) \ - * np.sin(np.pi * x[1]) * np.cos(np.pi * x[1]), - # - - np.exp(-8 * np.pi**2 * self.nu * t) * np.cos(2.0 * np.pi * (x[0] - t)) \ - * (np.cos(np.pi * x[1]))**2 - )) - ) - - w.sub(1).interpolate( lambda x: (4/17) * np.exp(-16 * np.pi**2 * self.nu * t) * \ - np.cos(4 * np.pi * (x[0]-t)) * np.cos(np.pi * x[1]) + 1 - ) - - + w.sub(0).interpolate( + lambda x: np.vstack( + ( + 1.0 + - np.exp(-8 * np.pi**2 * self.nu * t) + * np.sin(2.0 * np.pi * (x[0] - t)) + * np.sin(np.pi * x[1]) + * np.cos(np.pi * x[1]), + # + -np.exp(-8 * np.pi**2 * self.nu * t) + * np.cos(2.0 * np.pi * (x[0] - t)) + * (np.cos(np.pi * x[1])) ** 2, + ) + ) + ) + + w.sub(1).interpolate( + lambda x: (4 / 17) + * np.exp(-16 * np.pi**2 * self.nu * t) + * np.cos(4 * np.pi * (x[0] - t)) + * np.cos(np.pi * x[1]) + + 1 + ) me = self.dtype_u(self.init) self.convert_from_fenicsx_vector(input=w, output=me) @@ -459,7 +488,7 @@ def fix_residual(self, res): class fenicsx_NSE_periodic_mass(fenicsx_NSE_mass): - r""" + r""" Monolithic incompressible Navier–Stokes problem in 2D using FEniCSx. This class implements a finite element discretization of the incompressible @@ -475,7 +504,7 @@ class fenicsx_NSE_periodic_mass(fenicsx_NSE_mass): \quad \nabla \cdot \mathbf{u} = 0 Domain: - A square domain :math:`\Omega = [-0.5, 0.5]^2` discretized with triangular + A square domain :math:`\Omega = [-0.5, 0.5]^2` discretized with triangular elements. Spatial discretization: @@ -489,10 +518,10 @@ class fenicsx_NSE_periodic_mass(fenicsx_NSE_mass): u(x,y,t) = 1 - e^{-8\pi^{2}\nu t}\sin(2\pi(x-t))\sin(\pi y)\cos(\pi y) v(x,y,t) = - e^{-8\pi^{2}\nu t} \cos\(2\pi(x-t)) \cos^{2}(\pi y) p(x,y,t) = 1 + \frac{4}{17} e^{-16\pi^{2}\nu t} \cos\(4\pi(x-t)) \cos(\pi y) - - The solution is specifically constructed to investigate the impact of boundary - conditions on numerical accuracy. The source term g =(g_1, g_2) is derived from - the exact solution to ensure that it satisfies the Navier–Stokes equations. + + The solution is specifically constructed to investigate the impact of boundary + conditions on numerical accuracy. The source term g =(g_1, g_2) is derived from + the exact solution to ensure that it satisfies the Navier–Stokes equations. Boundary conditions: - periodic boundary conditions on the left and right boundaries @@ -501,7 +530,7 @@ class fenicsx_NSE_periodic_mass(fenicsx_NSE_mass): In this class, the same problem as "fenicsx_NSE_mass" is formulated with periodic boundary conditions on the left and right boundaries, while retaining the time- independent Dirichlet conditions on the top and bottom boundaries, in order to - isolate and compare the effects of time-dependent boundary data on the numerical + isolate and compare the effects of time-dependent boundary data on the numerical accuracy. @@ -550,49 +579,55 @@ def __init__(self, nelems=128, t0=0.0, family='CG', order=4, nu=0.1, comm=MPI.CO self.xmin = domain.geometry.x[:, 0].min() self.xmax = domain.geometry.x[:, 0].max() self.ymin = domain.geometry.x[:, 1].min() - self.ymax = domain.geometry.x[:, 1].max() - self.Lx = float(self.xmax - self.xmin) - self.Ly = float(self.ymax - self.ymin) + self.ymax = domain.geometry.x[:, 1].max() + self.Lx = float(self.xmax - self.xmin) + self.Ly = float(self.ymax - self.ymin) fdim = domain.topology.dim - 1 # locate facets on the right boundary for periodicity and create meshtags for them - facets_x = dfx.mesh.locate_entities_boundary(domain, fdim, self.periodic_boundary_x) + facets_x = dfx.mesh.locate_entities_boundary(domain, fdim, self.periodic_boundary_x) arg_sort_x = np.argsort(facets_x) mt_x = dfx.mesh.meshtags(domain, fdim, facets_x[arg_sort_x], np.full(len(facets_x), 2, dtype=np.int32)) - + # create multi-point constraints for periodicity and add the periodic constraints to tthe boundary conditions self.mpc = dolfinx_mpc.MultiPointConstraint(self.W) - self.mpc.create_periodic_constraint_topological(self.W.sub(0).sub(0), mt_x, 2, self.periodic_relation_x, self.bcs) - self.mpc.create_periodic_constraint_topological(self.W.sub(0).sub(1), mt_x, 2, self.periodic_relation_x, self.bcs) + self.mpc.create_periodic_constraint_topological( + self.W.sub(0).sub(0), mt_x, 2, self.periodic_relation_x, self.bcs + ) + self.mpc.create_periodic_constraint_topological( + self.W.sub(0).sub(1), mt_x, 2, self.periodic_relation_x, self.bcs + ) self.mpc.create_periodic_constraint_topological(self.W.sub(1), mt_x, 2, self.periodic_relation_x, self.bcs) self.mpc.finalize() - # reassemble the nonlinear problem + # reassemble the nonlinear problem self.factor = dfx.fem.Constant(domain, PETSc.ScalarType(0.0)) - + # redefine the solution function to be compatible with the multi-point constraints W_mpc = self.mpc.function_space self.w = dfx.fem.Function(W_mpc) u, p = ufl.split(self.w) rhs_u, rhs_p = ufl.split(self.rhs) - + # redefine the nonlinear form to be compatible with the multi-point constraints - F = (1/self.factor) * ufl.dot(u, self.v) * ufl.dx + F = (1 / self.factor) * ufl.dot(u, self.v) * ufl.dx F += ufl.dot(ufl.dot(u, ufl.nabla_grad(u)), self.v) * ufl.dx F += ufl.inner(self.nu * ufl.nabla_grad(u), ufl.nabla_grad(self.v)) * ufl.dx F -= ufl.dot(p, ufl.div(self.v)) * ufl.dx F -= ufl.dot(self.g, self.v) * ufl.dx F -= ufl.dot(ufl.div(u), self.q) * ufl.dx - F -= (1/self.factor) * ufl.dot(rhs_u, self.v) * ufl.dx - F -= (1/self.factor) * ufl.dot(rhs_p, self.q) * ufl.dx - + F -= (1 / self.factor) * ufl.dot(rhs_u, self.v) * ufl.dx + F -= (1 / self.factor) * ufl.dot(rhs_p, self.q) * ufl.dx + # compute the Jacobian of the nonlinear form for use in the Newton solver Jac = ufl.derivative(F, self.w, ufl.TrialFunction(self.W)) - + # assemble the nonlinear problem with the multi-point constraints and configure the Newton solver - self.problem = dolfinx_mpc.NonlinearProblem(F, self.w, mpc=self.mpc, bcs=self.bcs, J=Jac, petsc_options=self.petsc_options) - + self.problem = dolfinx_mpc.NonlinearProblem( + F, self.w, mpc=self.mpc, bcs=self.bcs, J=Jac, petsc_options=self.petsc_options + ) + def solve_system(self, rhs, factor, u0, t): r""" Dolfin's nonlinear solver for : @@ -629,7 +664,7 @@ def solve_system(self, rhs, factor, u0, t): self.w.x.scatter_forward() # self.convert_from_fenicsx_vector(input=self.w, output=u) - + return u def periodic_boundary_x(self, x): @@ -638,15 +673,15 @@ def periodic_boundary_x(self, x): Parameters ---------- x : array_like - Coordinates of the points to check for periodicity. + Coordinates of the points to check for periodicity. Returns ------- boolTrue if the point is on the right boundary (x = xmax), False otherwise. """ eps = 10 * np.finfo(x.dtype).resolution - return np.isclose(x[0], self.xmax, atol=eps) - + return np.isclose(x[0], self.xmax, atol=eps) + def periodic_relation_x(self, x): r""" Maps points from the right boundary to the left boundary for periodicity. @@ -664,5 +699,3 @@ def periodic_relation_x(self, x): out_x[1] = x[1] out_x[2] = x[2] return out_x - - \ No newline at end of file diff --git a/pySDC/projects/StroemungsRaum_FEniCSx/run_Navier_Stokes_equations_TaylGreen_monolithic_FEniCSx.py b/pySDC/projects/StroemungsRaum_FEniCSx/run_Navier_Stokes_equations_TaylGreen_monolithic_FEniCSx.py index 27f879f7a2..d784d914b5 100755 --- a/pySDC/projects/StroemungsRaum_FEniCSx/run_Navier_Stokes_equations_TaylGreen_monolithic_FEniCSx.py +++ b/pySDC/projects/StroemungsRaum_FEniCSx/run_Navier_Stokes_equations_TaylGreen_monolithic_FEniCSx.py @@ -5,9 +5,10 @@ from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from pySDC.projects.StroemungsRaum.sweepers.generic_implicit_mass import generic_implicit_mass from pySDC.projects.StroemungsRaum_FEniCSx.problem_classes.NavierStokes_2D_TaylGreen_monolithic_FEniCSx import ( - fenicsx_NSE_mass, - fenicsx_NSE_periodic_mass - ) + fenicsx_NSE_mass, + fenicsx_NSE_periodic_mass, +) + def setup(t0=0, periodic=False): """ @@ -45,27 +46,27 @@ def setup(t0=0, periodic=False): problem_params = dict() problem_params['nu'] = 0.2 - problem_params['t0'] = t0 + problem_params['t0'] = t0 problem_params['nelems'] = [64] problem_params['family'] = 'CG' problem_params['order'] = [2] # initialize controller parameters controller_params = dict() - controller_params['logger_level'] = 20 + controller_params['logger_level'] = 20 # Fill description dictionary for easy hierarchy creation description = dict() if periodic: - description['problem_class'] = fenicsx_NSE_periodic_mass + description['problem_class'] = fenicsx_NSE_periodic_mass else: - description['problem_class'] = fenicsx_NSE_mass + description['problem_class'] = fenicsx_NSE_mass description['sweeper_class'] = generic_implicit_mass description['problem_params'] = problem_params description['sweeper_params'] = sweeper_params description['level_params'] = level_params description['step_params'] = step_params - + return description, controller_params @@ -101,12 +102,13 @@ def run_simulation(description, controller_params, Tend): # call main function to get things done... uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - + return P, stats, uend + def run_postprocessing(P, uend, Tend): """ - Postprocess and store simulation results for visualization and analysis. + Compute L2 and relative L2 errors between numerical and exact solutions at final time. Args: P: Problem instance, @@ -116,7 +118,7 @@ def run_postprocessing(P, uend, Tend): Tend: float, Final simulation time. - Returns: + Returns: rel_error_u: float, Relative L2 error in velocity compared to the exact solution. rel_error_p: float, @@ -135,29 +137,29 @@ def run_postprocessing(P, uend, Tend): ue, pe = wx.split() ue_ = ue.collapse() pe_ = pe.collapse() - # + # error_u = np.sqrt(dfx.fem.assemble_scalar(dfx.fem.form(ufl.dot(un_ - ue_, un_ - ue_) * ufl.dx))) error_p = np.sqrt(dfx.fem.assemble_scalar(dfx.fem.form((pn_ - pe_) * (pn_ - pe_) * ufl.dx))) - - # + + # norm_u = np.sqrt(dfx.fem.assemble_scalar(dfx.fem.form(ufl.dot(ue_, ue_) * ufl.dx))) norm_p = np.sqrt(dfx.fem.assemble_scalar(dfx.fem.form(pe_ * pe_ * ufl.dx))) # rel_error_u = error_u / norm_u - rel_error_p = error_p / norm_p - + rel_error_p = error_p / norm_p + # print(f"L2 error in velocity : {error_u:.3e} (relative: {rel_error_u:.3e})") print(f"L2 error in pressure : {error_p:.3e} (relative: {rel_error_p:.3e})") return rel_error_u, rel_error_p - + if __name__ == "__main__": t0 = 0.0 - Tend =0.1 + Tend = 0.1 # run the setup to get description and controller parameters description, controller_params = setup(t0=t0, periodic=True) @@ -165,5 +167,5 @@ def run_postprocessing(P, uend, Tend): # run the simulation and get the problem, stats and final solution problem, stats, uend = run_simulation(description, controller_params, Tend) - # run postprocessing to save parameters and solutions for visualization - rel_error_u, rel_error_p = run_postprocessing(problem, uend, Tend) \ No newline at end of file + # run postprocessing to evaluate L2 errors + rel_error_u, rel_error_p = run_postprocessing(problem, uend, Tend) diff --git a/pySDC/projects/StroemungsRaum_FEniCSx/tests/test_Navier_Stokes_TaylGreen_monolithic_FEniCSx.py b/pySDC/projects/StroemungsRaum_FEniCSx/tests/test_Navier_Stokes_TaylGreen_monolithic_FEniCSx.py index a40ef04525..b228495937 100644 --- a/pySDC/projects/StroemungsRaum_FEniCSx/tests/test_Navier_Stokes_TaylGreen_monolithic_FEniCSx.py +++ b/pySDC/projects/StroemungsRaum_FEniCSx/tests/test_Navier_Stokes_TaylGreen_monolithic_FEniCSx.py @@ -5,13 +5,14 @@ import dolfinx.fem.petsc from pySDC.projects.StroemungsRaum_FEniCSx.problem_classes.NavierStokes_2D_TaylGreen_monolithic_FEniCSx import ( - fenicsx_NSE_mass, - fenicsx_NSE_periodic_mass, - ) + fenicsx_NSE_mass, + fenicsx_NSE_periodic_mass, +) + @pytest.mark.fenicsx def test_solve_system(): -# + # # Physical and numerical parameters for the test case nu = 0.02 t = 0.0 @@ -24,15 +25,15 @@ def test_solve_system(): # Evaluate the source term at the given time g = prob.source_term(t) - # get the exact solution at time t + # get the exact solution at time t wex = dfx.fem.Function(prob.W) wex.x.array[:] = prob.u_exact(t)[:] # Split the exact solution into velocity and pressure components u, p = ufl.split(wex) - + # assemble the right-hand side vector - rhs_f = ufl.dot(u, prob.v) * ufl.dx + rhs_f = ufl.dot(u, prob.v) * ufl.dx rhs_f += factor * ufl.dot(ufl.dot(u, ufl.nabla_grad(u)), prob.v) * ufl.dx rhs_f += factor * ufl.inner(prob.nu * ufl.nabla_grad(u), ufl.nabla_grad(prob.v)) * ufl.dx rhs_f -= factor * ufl.dot(p, ufl.div(prob.v)) * ufl.dx @@ -41,13 +42,12 @@ def test_solve_system(): b = dolfinx.fem.petsc.assemble_vector(dfx.fem.form(rhs_f)) b.assemble() - #b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) - + # b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # Solve the linear system for both non-periodic and periodic cases w = prob.solve_system(rhs=b, factor=factor, u0=prob.dtype_u(prob.init), t=t) w_periodic = prob_periodic.solve_system(rhs=b, factor=factor, u0=prob.dtype_u(prob.init), t=t) - + # get the exact solution at time t for error computation wx = prob.dtype_u(prob.init) wx[:] = wex.x.array[:] @@ -55,9 +55,10 @@ def test_solve_system(): # Compute the relative monolithic velocity–pressure error for both non-periodic and periodic cases rel_err = abs(wx - w) / abs(wx) rel_err_periodic = abs(wx - w_periodic) / abs(wx) - assert rel_err < 1e-5, (f"Relative monolithic velocity–pressure error {rel_err:.3e} exceeds tolerance") - assert rel_err_periodic < 1e-5, ( - f"Relative monolithic velocity–pressure error (periodic) {rel_err_periodic:.3e} exceeds tolerance") + assert rel_err < 1e-5, f"Relative monolithic velocity–pressure error {rel_err:.3e} exceeds tolerance" + assert ( + rel_err_periodic < 1e-5 + ), f"Relative monolithic velocity–pressure error (periodic) {rel_err_periodic:.3e} exceeds tolerance" @pytest.mark.fenicsx @@ -67,49 +68,65 @@ def test_eval_f(): nu = 0.2 t = 0.0 - # create the 2D monolithic Navier-Stokes problem + # create the 2D monolithic Navier-Stokes problem prob = fenicsx_NSE_mass(t0=0.0, nelems=64, family='CG', order=2, nu=nu) # get the exact solution at time t wex = dfx.fem.Function(prob.W) - wex.sub(0).interpolate(lambda x: - np.vstack((-1.0 *np.cos(x[0])*np.sin(x[1])*np.exp(-2.0*nu*t),np.sin(x[0])*np.cos(x[1])*np.exp(-2.0*nu*t)))) - wex.sub(1).interpolate( lambda x: -0.25*(np.cos(2.0*x[0]) + np.cos(2.0*x[1]))*np.exp(-4.0*nu*t)) + wex.sub(0).interpolate( + lambda x: np.vstack( + ( + -1.0 * np.cos(x[0]) * np.sin(x[1]) * np.exp(-2.0 * nu * t), + np.sin(x[0]) * np.cos(x[1]) * np.exp(-2.0 * nu * t), + ) + ) + ) + wex.sub(1).interpolate(lambda x: -0.25 * (np.cos(2.0 * x[0]) + np.cos(2.0 * x[1])) * np.exp(-4.0 * nu * t)) # evaluate the right-hand side vector using the problem's eval_f method f = prob.eval_f(w=wex.x.array[:], t=t) # compute the expected right-hand side vector by applying the mass matrix to the exact solution fw = dfx.fem.Function(prob.W) - fw.sub(0).interpolate(lambda x: np.vstack(( - 2.0*nu*np.cos(x[0])*np.sin(x[1])*np.exp(-2.0*nu*t) - + np.pi / 34 * np.exp(-16 * np.pi**2 * nu * t) * np.sin(4 * np.pi * (t - x[0])) - * np.cos(np.pi * x[1]) * (32 - 17*np.cos(np.pi*x[1])), - # - -2.0*nu*np.sin(x[0])*np.cos(x[1])*np.exp(-2.0*nu*t) - + 2 * np.pi**2 * nu * np.exp(-8 * np.pi**2 * nu * t) * np.cos(2 * np.pi * (t - x[0])) - - np.pi * np.exp(-16 * np.pi**2 * nu * t) * np.sin(np.pi * x[1]) * (2 * np.cos(np.pi \ - * x[1])**3 + 4/17 * np.cos(4 * np.pi * (t - x[0]))) - ))) - fw.sub(1).interpolate(lambda x: np.full(x.shape[1], 0.0, dtype=np.float64)) - + fw.sub(0).interpolate( + lambda x: np.vstack( + ( + 2.0 * nu * np.cos(x[0]) * np.sin(x[1]) * np.exp(-2.0 * nu * t) + + np.pi + / 34 + * np.exp(-16 * np.pi**2 * nu * t) + * np.sin(4 * np.pi * (t - x[0])) + * np.cos(np.pi * x[1]) + * (32 - 17 * np.cos(np.pi * x[1])), + # + -2.0 * nu * np.sin(x[0]) * np.cos(x[1]) * np.exp(-2.0 * nu * t) + + 2 * np.pi**2 * nu * np.exp(-8 * np.pi**2 * nu * t) * np.cos(2 * np.pi * (t - x[0])) + - np.pi + * np.exp(-16 * np.pi**2 * nu * t) + * np.sin(np.pi * x[1]) + * (2 * np.cos(np.pi * x[1]) ** 3 + 4 / 17 * np.cos(4 * np.pi * (t - x[0]))), + ) + ) + ) + fw.sub(1).interpolate(lambda x: np.full(x.shape[1], 0.0, dtype=np.float64)) + # apply the mass matrix to obtain the expected right-hand side vector Mfw = dfx.fem.Function(prob.W) prob.M.mult(fw.x.petsc_vec, Mfw.x.petsc_vec) # convert the expected right-hand side vector to a numpy array for comparison - fwx = prob.dtype_u(prob.init) + fwx = prob.dtype_u(prob.init) fwx[:] = Mfw.x.array[:] - + # apply boundary conditions to both vectors before comparison prob.fix_residual(f) prob.fix_residual(fwx) - + # compute the relative error between computed and expected right-hand sides rel_err = abs(f - fwx) / abs(fwx) print(f"Relative error in right-hand side evaluation: {rel_err:.3e}") assert rel_err < 1e-5, f"Relative error {rel_err} exceeds tolerance" - + @pytest.mark.fenicsx def test_problem_class(): @@ -117,7 +134,7 @@ def test_problem_class(): from pySDC.projects.StroemungsRaum_FEniCSx.run_Navier_Stokes_equations_TaylGreen_monolithic_FEniCSx import ( setup, run_simulation, - run_postprocessing + run_postprocessing, ) t0 = 0.0 @@ -132,7 +149,7 @@ def test_problem_class(): # run postprocessing to save parameters and solutions for visualization rel_error_u, rel_error_p = run_postprocessing(problem, uend, Tend) - + print(f"Relative L2 error in velocity: {rel_error_u:.3e}") print(f"Relative L2 error in pressure: {rel_error_p:.3e}") # assert that the relative errors in velocity and pressure are within acceptable limits @@ -148,7 +165,7 @@ def test_problem_class(): # run postprocessing to save parameters and solutions for visualization rel_error_u, rel_error_p = run_postprocessing(problem, uend, Tend) - + # assert that the relative errors in velocity and pressure are within acceptable limits assert rel_error_u < 1e-2, f"Error in velocity error (periodic bc) {rel_error_u} exceeds tolerance" - assert rel_error_p < 1e-2, f"Error in pressure error (periodic bc) {rel_error_p} exceeds tolerance" \ No newline at end of file + assert rel_error_p < 1e-2, f"Error in pressure error (periodic bc) {rel_error_p} exceeds tolerance"