Coverage for pySDC / projects / StroemungsRaum / problem_classes / ConvectionDiffusion_2D_FEniCS.py: 100%
58 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-27 07:06 +0000
1import dolfin as df
2from pySDC.core.problem import Problem
3from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh
6class fenics_ConvDiff2D_mass(Problem):
7 r"""
8 Example implementing a forced two-dimensional convection-diffusion equation using Dirichlet
9 boundary conditions. The problem considered is a benchmark test with a rotating Gaussian profile.
11 The equation we are solving is the two-dimensional convection-diffusion equation:
13 .. math::
14 \frac{\partial u}{\partial t} = - U \cdot \nabla u + \nu \Delta u + f
16 where:
17 .. math:`u(x, y, t)` is the scalar field we are solving for (e.g., concentration or temperature).
18 .. math:`U = (u, v)` is the velocity field of the flow.
19 .. math:`\nu` is the kinematic viscosity, which quantifies the diffusion rate.
20 .. math:`f(x, y, t)` is the source term, representing external forcing or generation of .. math:`u`.
22 The computational domain for this problem is:
24 .. math::
25 x \in \Omega := [-0.5, 0.5] \times [-0.5, 0.5]
27 Dirichlet boundary conditions are applied, meaning that the value of .. math:`u` is specified on the
28 boundary of the domain. In this benchmark example, the forcing term .. math:`f` is:
30 .. math::
31 f(x,y,t) = 0
33 This implies there are no additional sources or sinks affecting the field .. math:`u`, simplifying the problem to just the
34 effects of convection and diffusion. The analytical solution for the scalar field .. math:`u is given by:
36 .. math::
37 u(x,y,t) = \frac{\sigma^2}{\sigma^2 + 4 \nu t} \exp\left( - \frac{(\hat{x}-x_0)^2 + (\hat{y}-y_0)^2}{\sigma^2 + 4 \nu t} \right)
39 where:
40 .. math:`\sigma` is the initial standard deviation of the Gaussian.
41 .. math:`\omega` is the angular velocity of the rotation (in this example, :math:`\omega = 4`).
42 .. math:`(\hat{x}, \hat{y})` are the rotated coordinates defined by:
43 .. math::
44 \hat{x} = \cos(4 t) x + \sin(4 t) y \\
45 \hat{y} = -\sin(4 t) x + \cos(4 t) y
46 .. math:`(x_0, y_0)` is the center of the Gaussian, given as .. math:`(-0.25, 0.0)`.
49 The velocity field .. math:`U` for the rotating Gaussian is defined as:
51 .. math::
52 U = (u,v) = (-4*y, 4*x)
54 This represents a counter-clockwise rotation around the origin with angular velocity .. math:`\omega`.
55 The flow causes the scalar field .. math:`u` to be advected in a circular pattern, simulating the rotation of the Gaussian.
58 Parameters
59 ----------
60 c_nvars : int, optional
61 Spatial resolution, i.e., numbers of degrees of freedom in space.
62 t0 : float, optional
63 Starting time.
64 family : str, optional
65 Indicates the family of elements used to create the function space
66 for the trail and test functions. The default is ``'CG'``, which are the class
67 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [2]_.
68 order : int, optional
69 Defines the order of the elements in the function space.
70 sigma : float, optional
71 Coefficient associated with the mass term or reaction term in the equation.
72 nu : float, optional
73 Diffusion coefficient :math:`\nu`.
75 Attributes
76 ----------
77 V : FunctionSpace
78 Defines the function space of the trial and test functions.
79 M : scalar, vector, matrix or higher rank tensor
80 Denotes the expression :math:`\int_\Omega u_t v\,dx`.
81 K : scalar, vector, matrix or higher rank tensor
82 Denotes the expression :math:`- \nu \int_\Omega \nabla u \nabla v\,dx`.
83 g : Expression
84 The forcing term :math:`f` in the convection-diffusion equations.
85 bc : DirichletBC
86 Denotes the Dirichlet boundary conditions.
88 References
89 ----------
90 .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
91 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
92 .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
93 Wells and others. Springer (2012).
94 """
96 dtype_u = fenics_mesh
97 dtype_f = rhs_fenics_mesh
99 def __init__(self, c_nvars=64, t0=0.0, family='CG', order=2, nu=0.01, sigma=0.05):
101 # define the Dirichlet boundary
102 def Boundary(x, on_boundary):
103 return on_boundary
105 # set solver and form parameters
106 df.parameters["form_compiler"]["optimize"] = True
107 df.parameters["form_compiler"]["cpp_optimize"] = True
108 df.parameters['allow_extrapolation'] = True
110 # set mesh and refinement (for multilevel)
111 mesh = df.RectangleMesh(df.Point(-0.5, -0.5), df.Point(0.5, 0.5), c_nvars, c_nvars)
113 # define function space for future reference
114 self.V = df.FunctionSpace(mesh, family, order)
115 tmp = df.Function(self.V)
116 self.logger.debug('DoFs on this level:', len(tmp.vector()[:]))
118 # define velocity
119 self.VC = df.VectorFunctionSpace(mesh, family, order)
120 self.U = df.interpolate(df.Expression(('-4*x[1]', '4*x[0]'), degree=order), self.VC)
122 super().__init__(self.V)
123 self._makeAttributeAndRegister(
124 'c_nvars', 't0', 'family', 'order', 'nu', 'sigma', localVars=locals(), readOnly=True
125 )
127 # define trial and test functions
128 u = df.TrialFunction(self.V)
129 v = df.TestFunction(self.V)
131 # mass, stiffness and convection terms
132 a_K = -1.0 * df.inner(df.nabla_grad(u), self.nu * df.nabla_grad(v)) * df.dx
133 a_M = u * v * df.dx
134 a_C = -1.0 * df.inner(df.dot(self.U, df.nabla_grad(u)), v) * df.dx
136 # assemble the forms
137 self.M = df.assemble(a_M)
138 self.K = df.assemble(a_K)
139 self.C = df.assemble(a_C)
141 # set exact solution for boundary conditions
142 self.u_D = df.Expression(
143 'pow(s,2)/(pow(s,2)+4*nu*t)*exp(-(pow(((cos(4*t)*x[0]+sin(4*t)*x[1])-x0),2)\
144 +pow(((-sin(4*t)*x[0]+cos(4*t)*x[1])-y0),2))/(pow(s,2)+4*nu*t))',
145 s=self.sigma,
146 nu=self.nu,
147 x0=-0.25,
148 y0=0.0,
149 t=t0,
150 degree=self.order,
151 )
153 # set boundary conditions
154 self.bc = df.DirichletBC(self.V, self.u_D, Boundary)
155 self.bc_hom = df.DirichletBC(self.V, df.Constant(0), Boundary)
156 self.fix_bc_for_residual = True
158 # set forcing term as expression
159 self.g = df.Expression('0.0', t=self.t0, degree=self.order)
161 def solve_system(self, rhs, factor, u0, t):
162 r"""
163 Dolfin's linear solver for :math:`(M - factor A) \vec{u} = \vec{rhs}`.
165 Parameters
166 ----------
167 rhs : dtype_f
168 Right-hand side for the nonlinear system.
169 factor : float
170 Abbrev. for the node-to-node stepsize (or any other factor required).
171 u0 : dtype_u
172 Initial guess for the iterative solver (not used here so far).
173 t : float
174 Current time.
176 Returns
177 -------
178 u : dtype_u
179 Solution.
180 """
181 self.u_D.t = t
183 u = self.dtype_u(u0)
184 T = self.M - factor * self.K
186 self.bc.apply(T, rhs.values.vector())
187 df.solve(T, u.values.vector(), rhs.values.vector())
189 return u
191 def eval_f(self, u, t):
192 """
193 Routine to evaluate both parts of the right-hand side of the problem.
195 Parameters
196 ----------
197 u : dtype_u
198 Current values of the numerical solution.
199 t : float
200 Current time at which the numerical solution is computed.
202 Returns
203 -------
204 f : dtype_f
205 The right-hand side divided into two parts.
206 """
207 self.g.t = t
209 f = self.dtype_f(self.V)
210 self.K.mult(u.values.vector(), f.impl.values.vector())
211 self.C.mult(u.values.vector(), f.expl.values.vector())
213 f.expl += self.apply_mass_matrix(self.dtype_u(df.interpolate(self.g, self.V)))
215 return f
217 def u_exact(self, t):
218 r"""
219 Routine to compute the exact solution at time :math:`t`.
221 Parameters
222 ----------
223 t : float
224 Time of the exact solution.
226 Returns
227 -------
228 me : dtype_u
229 Exact solution.
230 """
232 u0 = df.Expression(
233 'pow(s,2)/(pow(s,2)+4*nu*t)*exp(-(pow(((cos(4*t)*x[0]+sin(4*t)*x[1])-x0),2)\
234 +pow(((-sin(4*t)*x[0]+cos(4*t)*x[1])-y0),2))/(pow(s,2)+4*nu*t))',
235 s=self.sigma,
236 nu=self.nu,
237 x0=-0.25,
238 y0=0.0,
239 t=t,
240 degree=self.order,
241 )
243 me = self.dtype_u(df.interpolate(u0, self.V))
245 return me
247 def apply_mass_matrix(self, u):
248 r"""
249 Routine to apply mass matrix.
251 Parameters
252 ----------
253 u : dtype_u
254 Current values of the numerical solution.
256 Returns
257 -------
258 me : dtype_u
259 The product :math:`M \vec{u}`.
260 """
262 me = self.dtype_u(self.V)
263 self.M.mult(u.values.vector(), me.values.vector())
265 return me
267 def fix_residual(self, res):
268 """
269 Applies homogeneous Dirichlet boundary conditions to the residual
271 Parameters
272 ----------
273 res : dtype_u
274 Residual
275 """
276 self.bc_hom.apply(res.values.vector())
277 return None