|
| 1 | +import dolfin as df |
| 2 | +from pySDC.core.problem import Problem |
| 3 | +from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh |
| 4 | + |
| 5 | + |
| 6 | +class 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. |
| 10 | +
|
| 11 | + The equation we are solving is the two-dimensional convection-diffusion equation: |
| 12 | +
|
| 13 | + .. math:: |
| 14 | + \frac{\partial u}{\partial t} = - U \cdot \nabla u + \nu \Delta u + f |
| 15 | +
|
| 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`. |
| 21 | +
|
| 22 | + The computational domain for this problem is: |
| 23 | +
|
| 24 | + .. math:: |
| 25 | + x \in \Omega := [-0.5, 0.5] \times [-0.5, 0.5] |
| 26 | +
|
| 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: |
| 29 | +
|
| 30 | + .. math:: |
| 31 | + f(x,y,t) = 0 |
| 32 | +
|
| 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: |
| 35 | +
|
| 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) |
| 38 | +
|
| 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)`. |
| 47 | +
|
| 48 | +
|
| 49 | + The velocity field .. math:`U` for the rotating Gaussian is defined as: |
| 50 | +
|
| 51 | + .. math:: |
| 52 | + U = (u,v) = (-4*y, 4*x) |
| 53 | +
|
| 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. |
| 56 | +
|
| 57 | +
|
| 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`. |
| 74 | +
|
| 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. |
| 87 | +
|
| 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 | + """ |
| 95 | + |
| 96 | + dtype_u = fenics_mesh |
| 97 | + dtype_f = rhs_fenics_mesh |
| 98 | + |
| 99 | + def __init__(self, c_nvars=64, t0=0.0, family='CG', order=2, nu=0.01, sigma=0.05): |
| 100 | + |
| 101 | + # define the Dirichlet boundary |
| 102 | + def Boundary(x, on_boundary): |
| 103 | + return on_boundary |
| 104 | + |
| 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 |
| 109 | + |
| 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) |
| 112 | + |
| 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()[:])) |
| 117 | + |
| 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) |
| 121 | + |
| 122 | + super().__init__(self.V) |
| 123 | + self._makeAttributeAndRegister( |
| 124 | + 'c_nvars', 't0', 'family', 'order', 'nu', 'sigma', localVars=locals(), readOnly=True |
| 125 | + ) |
| 126 | + |
| 127 | + # define trial and test functions |
| 128 | + u = df.TrialFunction(self.V) |
| 129 | + v = df.TestFunction(self.V) |
| 130 | + |
| 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 |
| 135 | + |
| 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) |
| 140 | + |
| 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 | + ) |
| 152 | + |
| 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 |
| 157 | + |
| 158 | + # set forcing term as expression |
| 159 | + self.g = df.Expression('0.0', t=self.t0, degree=self.order) |
| 160 | + |
| 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}`. |
| 164 | +
|
| 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. |
| 175 | +
|
| 176 | + Returns |
| 177 | + ------- |
| 178 | + u : dtype_u |
| 179 | + Solution. |
| 180 | + """ |
| 181 | + self.u_D.t = t |
| 182 | + |
| 183 | + u = self.dtype_u(u0) |
| 184 | + T = self.M - factor * self.K |
| 185 | + |
| 186 | + self.bc.apply(T, rhs.values.vector()) |
| 187 | + df.solve(T, u.values.vector(), rhs.values.vector()) |
| 188 | + |
| 189 | + return u |
| 190 | + |
| 191 | + def eval_f(self, u, t): |
| 192 | + """ |
| 193 | + Routine to evaluate both parts of the right-hand side of the problem. |
| 194 | +
|
| 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. |
| 201 | +
|
| 202 | + Returns |
| 203 | + ------- |
| 204 | + f : dtype_f |
| 205 | + The right-hand side divided into two parts. |
| 206 | + """ |
| 207 | + self.g.t = t |
| 208 | + |
| 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()) |
| 212 | + |
| 213 | + f.expl += self.apply_mass_matrix(self.dtype_u(df.interpolate(self.g, self.V))) |
| 214 | + |
| 215 | + return f |
| 216 | + |
| 217 | + def u_exact(self, t): |
| 218 | + r""" |
| 219 | + Routine to compute the exact solution at time :math:`t`. |
| 220 | +
|
| 221 | + Parameters |
| 222 | + ---------- |
| 223 | + t : float |
| 224 | + Time of the exact solution. |
| 225 | +
|
| 226 | + Returns |
| 227 | + ------- |
| 228 | + me : dtype_u |
| 229 | + Exact solution. |
| 230 | + """ |
| 231 | + |
| 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 | + ) |
| 242 | + |
| 243 | + me = self.dtype_u(df.interpolate(u0, self.V)) |
| 244 | + |
| 245 | + return me |
| 246 | + |
| 247 | + def apply_mass_matrix(self, u): |
| 248 | + r""" |
| 249 | + Routine to apply mass matrix. |
| 250 | +
|
| 251 | + Parameters |
| 252 | + ---------- |
| 253 | + u : dtype_u |
| 254 | + Current values of the numerical solution. |
| 255 | +
|
| 256 | + Returns |
| 257 | + ------- |
| 258 | + me : dtype_u |
| 259 | + The product :math:`M \vec{u}`. |
| 260 | + """ |
| 261 | + |
| 262 | + me = self.dtype_u(self.V) |
| 263 | + self.M.mult(u.values.vector(), me.values.vector()) |
| 264 | + |
| 265 | + return me |
| 266 | + |
| 267 | + def fix_residual(self, res): |
| 268 | + """ |
| 269 | + Applies homogeneous Dirichlet boundary conditions to the residual |
| 270 | +
|
| 271 | + Parameters |
| 272 | + ---------- |
| 273 | + res : dtype_u |
| 274 | + Residual |
| 275 | + """ |
| 276 | + self.bc_hom.apply(res.values.vector()) |
| 277 | + return None |
0 commit comments