Coverage for pySDC/projects/StroemungsRaum/problem_classes/ConvectionDiffusion_2D_FEniCS.py: 100%
59 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +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 use_restricted_boundary_values(self, u, t):
162 """
163 Refuse multilevel runs: the solves impose the time-dependent boundary data, which at the coarse nodes is not
164 the restricted fine boundary values, so MLSDC and PFASST would not converge to the collocation solution. See
165 :meth:`pySDC.core.problem.Problem.use_restricted_boundary_values`.
166 """
167 raise NotImplementedError(f'{type(self).__name__} cannot impose restricted boundary values on coarse levels')
169 def solve_system(self, rhs, factor, u0, t):
170 r"""
171 Dolfin's linear solver for :math:`(M - factor A) \vec{u} = \vec{rhs}`.
173 Parameters
174 ----------
175 rhs : dtype_f
176 Right-hand side for the nonlinear system.
177 factor : float
178 Abbrev. for the node-to-node stepsize (or any other factor required).
179 u0 : dtype_u
180 Initial guess for the iterative solver (not used here so far).
181 t : float
182 Current time.
184 Returns
185 -------
186 u : dtype_u
187 Solution.
188 """
189 self.u_D.t = t
191 u = self.dtype_u(u0)
192 T = self.M - factor * self.K
194 self.bc.apply(T, rhs.values.vector())
195 df.solve(T, u.values.vector(), rhs.values.vector())
197 return u
199 def eval_f(self, u, t):
200 """
201 Routine to evaluate both parts of the right-hand side of the problem.
203 Parameters
204 ----------
205 u : dtype_u
206 Current values of the numerical solution.
207 t : float
208 Current time at which the numerical solution is computed.
210 Returns
211 -------
212 f : dtype_f
213 The right-hand side divided into two parts.
214 """
215 self.g.t = t
217 f = self.dtype_f(self.V)
218 self.K.mult(u.values.vector(), f.impl.values.vector())
219 self.C.mult(u.values.vector(), f.expl.values.vector())
221 f.expl += self.apply_mass_matrix(self.dtype_u(df.interpolate(self.g, self.V)))
223 return f
225 def u_exact(self, t):
226 r"""
227 Routine to compute the exact solution at time :math:`t`.
229 Parameters
230 ----------
231 t : float
232 Time of the exact solution.
234 Returns
235 -------
236 me : dtype_u
237 Exact solution.
238 """
240 u0 = df.Expression(
241 'pow(s,2)/(pow(s,2)+4*nu*t)*exp(-(pow(((cos(4*t)*x[0]+sin(4*t)*x[1])-x0),2)\
242 +pow(((-sin(4*t)*x[0]+cos(4*t)*x[1])-y0),2))/(pow(s,2)+4*nu*t))',
243 s=self.sigma,
244 nu=self.nu,
245 x0=-0.25,
246 y0=0.0,
247 t=t,
248 degree=self.order,
249 )
251 me = self.dtype_u(df.interpolate(u0, self.V))
253 return me
255 def apply_mass_matrix(self, u):
256 r"""
257 Routine to apply mass matrix.
259 Parameters
260 ----------
261 u : dtype_u
262 Current values of the numerical solution.
264 Returns
265 -------
266 me : dtype_u
267 The product :math:`M \vec{u}`.
268 """
270 me = self.dtype_u(self.V)
271 self.M.mult(u.values.vector(), me.values.vector())
273 return me
275 def fix_residual(self, res):
276 """
277 Applies homogeneous Dirichlet boundary conditions to the residual
279 Parameters
280 ----------
281 res : dtype_u
282 Residual
283 """
284 self.bc_hom.apply(res.values.vector())
285 return None