|
| 1 | +Steady Boussinesq problem with integral constraints |
| 2 | +=================================================== |
| 3 | + |
| 4 | +.. rst-class:: emphasis |
| 5 | + |
| 6 | + This demo demonstrates how integral constraints can be imposed on a field |
| 7 | + that appears nonlinearly in the governing equations. |
| 8 | + |
| 9 | + The demo was contributed by `Aaron Baier-Reinio |
| 10 | + <mailto:[email protected]>`__ and ` Kars Knook |
| 11 | + |
| 12 | + |
| 13 | +We consider a nondimensionalised steady Boussinesq problem. |
| 14 | +The domain is :math:`\Omega \subset \mathbb{R}^2` |
| 15 | +with boundary :math:`\Gamma`. The boundary value problem is: |
| 16 | + |
| 17 | +.. math:: |
| 18 | +
|
| 19 | + -\nabla \cdot (\mu(T) \epsilon (u)) + \nabla p - f T &= 0 \quad \textrm{in}\ \Omega, |
| 20 | +
|
| 21 | + \nabla \cdot u &= 0 \quad \textrm{in}\ \Omega, |
| 22 | +
|
| 23 | + u \cdot \nabla T - \nabla^2 T &= 0 \quad \textrm{in}\ \Omega, |
| 24 | +
|
| 25 | + u &= 0 \quad \textrm{on}\ \Gamma, |
| 26 | +
|
| 27 | + \nabla T \cdot n &= g \quad \textrm{on}\ \Gamma. |
| 28 | +
|
| 29 | +The unknowns are the :math:`\mathbb{R}^2`-valued velocity :math:`u`, |
| 30 | +scalar pressure :math:`p` and temperature :math:`T`. |
| 31 | +The viscosity :math:`\mu(T)` is assumed to be a known function of :math:`T`. |
| 32 | +Moreover :math:`\epsilon (u)` denotes the symmetric gradient of :math:`u` |
| 33 | +and :math:`f = (0, -1)^T` the acceleration due to gravity. |
| 34 | +Inhomogeneous Neumann boundary conditions are imposed on :math:`\nabla T \cdot n` through |
| 35 | +given data :math:`g` which must satisfy a compatibility condition |
| 36 | + |
| 37 | +.. math:: |
| 38 | +
|
| 39 | + \int_{\Gamma} g \ {\rm d} s = 0. |
| 40 | +
|
| 41 | +Evidently the pressure :math:`p` is only determined up to a constant, since it only appears in |
| 42 | +the problem through its gradient. This choice of constant is arbitrary and does not affect the model. |
| 43 | +The situation regarding the temperature :math:`T` is, however, more subtle. |
| 44 | +For the sake of discussion let us first assume that :math:`\mu(T) = \mu_0` is a constant that does |
| 45 | +not depend on :math:`T`. It is then clear that, just like the pressure, the temperature :math:`T` |
| 46 | +is undetermined up to a constant. We shall pin this down by enforcing that the mean of :math:`T` |
| 47 | +is a user-supplied constant :math:`T_0`, |
| 48 | + |
| 49 | +.. math:: |
| 50 | +
|
| 51 | + \int_{\Omega} (T - T_0) \ {\rm d} x = 0. |
| 52 | +
|
| 53 | +The Boussinesq approximation assumes that density varies linearly with temperature. |
| 54 | +Hence, this constraint can be viewed as an indirect imposition on the total mass of fluid in :math:`\Omega`. |
| 55 | + |
| 56 | +Now suppose that :math:`\mu(T)` does depend on :math:`T`. |
| 57 | +For simplicity we use a basic power law :math:`\mu(T) = 10^{-4} \cdot T^{1/2}` |
| 58 | +but emphasise that the precise functional form of :math:`\mu(T)` is unimportant to this demo. |
| 59 | +We must still impose the integral constraint on :math:`T` to obtain a unique solution, |
| 60 | +but the value of :math:`T_0` now affects the model in a non-trivial way since :math:`\mu(T)` and |
| 61 | +:math:`T` are coupled (c.f. the figures at the bottom of the demo). |
| 62 | +**In particular, this is not a "trivial" situation like the incompressible |
| 63 | +Stokes problem where the discretized Jacobians have a nullspace corresponding to the constant pressures. |
| 64 | +Instead, we have an integral constraint on** :math:`T` **even though |
| 65 | +the discretized Jacobians do not have a nullspace corresponding to the constant temperatures.** |
| 66 | + |
| 67 | +We build the mesh using :doc:`netgen <netgen_mesh.py>`, choosing a trapezoidal geometry |
| 68 | +to prevent hydrostatic equilibrium and allow for a non-trivial velocity solution.:: |
| 69 | + |
| 70 | + from firedrake import * |
| 71 | + import netgen.occ as ngocc |
| 72 | + |
| 73 | + wp = ngocc.WorkPlane() |
| 74 | + wp.MoveTo(0, 0) |
| 75 | + wp.LineTo(2, 0) |
| 76 | + wp.LineTo(2, 1) |
| 77 | + wp.LineTo(0, 2) |
| 78 | + wp.LineTo(0, 0) |
| 79 | + |
| 80 | + shape = wp.Face() |
| 81 | + shape.edges.Min(ngocc.X).name = "left" |
| 82 | + shape.edges.Max(ngocc.X).name = "right" |
| 83 | + |
| 84 | + ngmesh = ngocc.OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.1) |
| 85 | + |
| 86 | + left_id = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(dim=1)) if name == "left"] |
| 87 | + right_id = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(dim=1)) if name == "right"] |
| 88 | + |
| 89 | + mesh = Mesh(ngmesh) |
| 90 | + x, y = SpatialCoordinate(mesh) |
| 91 | + |
| 92 | +Next we set up the discrete function spaces. |
| 93 | +We use lowest-order Taylor--Hood elements for the velocity and pressure, |
| 94 | +and continuous piecewise-linear elements for the temperature. |
| 95 | +We introduce a Lagrange multiplier to enforce the integral constraint on :math:`T`:: |
| 96 | + |
| 97 | + U = VectorFunctionSpace(mesh, "CG", degree=2) |
| 98 | + V = FunctionSpace(mesh, "CG", degree=1) |
| 99 | + W = FunctionSpace(mesh, "CG", degree=1) |
| 100 | + R = FunctionSpace(mesh, "R", 0) |
| 101 | + |
| 102 | + Z = U * V * W * R |
| 103 | + |
| 104 | +The trial and test functions are::: |
| 105 | + |
| 106 | + z = Function(Z) |
| 107 | + (u, p, T_aux, l) = split(z) |
| 108 | + (v, q, w, s) = split(TestFunction(Z)) |
| 109 | + |
| 110 | + T = T_aux + l |
| 111 | + |
| 112 | +The test Lagrange multiplier :code:`s` will allow us to impose the integral constraint on the temperature. |
| 113 | +We use the trial Lagrange multiplier :code:`l` by decomposing the discretized temperature field :code:`T` |
| 114 | +as :code:`T = T_aux + l` where :code:`T_aux` is the trial function from :code:`W`. |
| 115 | +The value of :code:`l` will then be determined by the integral constraint on :code:`T`. |
| 116 | + |
| 117 | +The remaining problem data to be specified is the Neumann data, |
| 118 | +viscosity, acceleration due to gravity and :math:`T_0`. |
| 119 | +For the Neumann data we choose a parabolic profile on the left and right edges, |
| 120 | +and zero data on the top and bottom.:: |
| 121 | + |
| 122 | + g_left = y*(y-2) # Neumann data on the left |
| 123 | + g_right = -8*y*(y-1) # Neumann data on the right |
| 124 | + mu = (1e-4) * T ** (1/2) # Viscosity |
| 125 | + f = as_vector([0, -1]) # Acceleration due to gravity |
| 126 | + T0 = Constant(1) # Mean of the temperature |
| 127 | + |
| 128 | +The nonlinear form for the problem is::: |
| 129 | + |
| 130 | + F = (mu * inner(sym(grad(u)), sym(grad(v))) * dx # Viscous terms |
| 131 | + - inner(p, div(v)) * dx # Pressure gradient |
| 132 | + - inner(f, T*v) * dx # Buoyancy term |
| 133 | + - inner(div(u), q) * dx # Incompressibility constraint |
| 134 | + + inner(dot(u, grad(T)), w) * dx # Temperature advection term |
| 135 | + + inner(grad(T), grad(w)) * dx # Temperature diffusion term |
| 136 | + - inner(g_left, w) * ds(tuple(left_id)) # Weakly imposed Neumann BC terms |
| 137 | + - inner(g_right, w) * ds(tuple(right_id)) # Weakly imposed Neumann BC terms |
| 138 | + + inner(T - T0, s) * dx # Integral constraint on T |
| 139 | + ) |
| 140 | + |
| 141 | +and the (strongly enforced) Dirichlet boundary conditions on :math:`u` are enforced by:: |
| 142 | + |
| 143 | + bc_u = DirichletBC(Z.sub(0), 0, "on_boundary") |
| 144 | + |
| 145 | +At this point we could form and solve a :class:`~.NonlinearVariationalProblem` |
| 146 | +using :code:`F` and :code:`bc_u`. However, the resultant problem has a nullspace of |
| 147 | +dimension 2, corresponding to (i) shifting :math:`p` by a constant :math:`C_1` |
| 148 | +and (ii) shifting :math:`l` by a constant :math:`C_2` while simultaneuosly shifting |
| 149 | +:math:`T_{\textrm{aux}}` by :math:`-C_2`. |
| 150 | + |
| 151 | +One way of dealing with nullspaces in Firedrake is to pass a :code:`nullspace` and |
| 152 | +:code:`transpose_nullspace` to :class:`~.NonlinearVariationalSolver`. However, sometimes |
| 153 | +this approach may not be practical. First, for nonlinear problems with Jacobians that |
| 154 | +are not symmetric, it may not obvious what the :code:`transpose_nullspace` is. A second |
| 155 | +reason is that, when using customised PETSc linear solvers, it may be desirable |
| 156 | +to directly eliminate the nullspace from the assembled Jacobian matrix, since one |
| 157 | +cannot always be sure that the linear solver at hand is correctly utilising the provided |
| 158 | +:code:`nullspace` and :code:`transpose_nullspace`. |
| 159 | + |
| 160 | +To directly eliminate the nullspace we introduce a class :code:`FixAtPointBC` which |
| 161 | +implements a boundary condition that fixes a field at a single point.:: |
| 162 | + |
| 163 | + import firedrake.utils as firedrake_utils |
| 164 | + |
| 165 | + class FixAtPointBC(firedrake.DirichletBC): |
| 166 | + r'''A special BC object for pinning a function at a point. |
| 167 | + |
| 168 | + :arg V: the :class:`.FunctionSpace` on which the boundary condition should be applied. |
| 169 | + :arg g: the boundary condition value. |
| 170 | + :arg bc_point: the point at which to pin the function. |
| 171 | + The location of the finite element DOF nearest to bc_point is actually used. |
| 172 | + ''' |
| 173 | + def __init__(self, V, g, bc_point): |
| 174 | + super(FixAtPointBC, self).__init__(V, g, bc_point) |
| 175 | + if isinstance(bc_point, tuple): |
| 176 | + bc_point = as_vector(bc_point) |
| 177 | + self.bc_point = bc_point |
| 178 | + |
| 179 | + @firedrake_utils.cached_property |
| 180 | + def nodes(self): |
| 181 | + V = self.function_space() |
| 182 | + x = firedrake.SpatialCoordinate(V.mesh()) |
| 183 | + xdist = x - self.bc_point |
| 184 | + |
| 185 | + test = firedrake.TestFunction(V) |
| 186 | + trial = firedrake.TrialFunction(V) |
| 187 | + xphi = firedrake.assemble(ufl.inner(xdist * test, xdist * trial) * ufl.dx, diagonal=True) |
| 188 | + phi = firedrake.assemble(ufl.inner(test, trial) * ufl.dx, diagonal=True) |
| 189 | + with xphi.dat.vec as xu, phi.dat.vec as u: |
| 190 | + xu.pointwiseDivide(xu, u) |
| 191 | + min_index, min_value = xu.min() # Find the index of the DOF closest to bc_point |
| 192 | + |
| 193 | + nodes = V.dof_dset.lgmap.applyInverse([min_index]) |
| 194 | + nodes = nodes[nodes >= 0] |
| 195 | + return nodes |
| 196 | + |
| 197 | +We use this to fix the pressure and auxiliary temperature at the origin::: |
| 198 | + |
| 199 | + aux_bcs = [FixAtPointBC(Z.sub(1), 0, as_vector([0, 0])), |
| 200 | + FixAtPointBC(Z.sub(2), 0, as_vector([0, 0]))] |
| 201 | + |
| 202 | +:code:`FixAtPointBC` takes three arguments: the function space to fix, the value with which it |
| 203 | +will be fixed, and the location at which to fix. Generally :code:`FixAtPointBC` will not fix |
| 204 | +the function at exactly the supplied point; it will fix it at the finite element DOF closest |
| 205 | +to that point. By default CG elements have DOFs on all mesh vertices, so if the supplied |
| 206 | +point if a mesh vertex then CG fields will be fixed at exactly the supplied point. |
| 207 | + |
| 208 | +.. warning:: |
| 209 | + |
| 210 | + A :code:`FixAtPointBC` does more than just fix the corresponding trial function |
| 211 | + at the chosen DOF. It also ensures that the corresponding test function |
| 212 | + (which is equal to one at that DOF and zero at all others) |
| 213 | + will no longer be used in the discretized variational problem. |
| 214 | + One must be sure that it is mathematically acceptable to do this. |
| 215 | + |
| 216 | + In the present setting this is acceptable owing to the homogeneous Dirichlet |
| 217 | + boundary conditions on :math:`u` and compatibility condition :math:`\int_{\Gamma} g \ {\rm d} s = 0` |
| 218 | + on the Neumann data. The former ensures that the rows in the discretized Jacobian |
| 219 | + corresponding to the incompressibility constraint are linearly dependent |
| 220 | + (if there are :math:`M` rows, only :math:`M-1` of them are linearly independent, since |
| 221 | + the boundary conditions on :math:`u` ensure that |
| 222 | + :math:`\int_{\Omega} \nabla \cdot u {\rm d} x = 0` automatically). |
| 223 | + Similarily the rows in the Jacobian corresponding to the temperature advection-diffusion |
| 224 | + equation are linearly independent (again, if there are :math:`M` rows, |
| 225 | + only :math:`M-1` of them are linearly independent). |
| 226 | + The effect of :code:`FixAtPointBC` will be to remove one of the rows corresponding |
| 227 | + to the incompressibility constraint and one corresponding to the temperature advection-diffusion |
| 228 | + equation. Which row ends up getting removed is determined by the location of :code:`bc_point`, |
| 229 | + but in the present setting removing a given row is mathematically equivalent to removing any one of the others. |
| 230 | + |
| 231 | + One could envisage a more complicated scenario than the one in this demo, wherein the Neumann data |
| 232 | + depends nonlinearly on some other problem unknowns, and it only satisfies the compatibility condition |
| 233 | + approximately (e.g. up to some discretization error). |
| 234 | + In this case one would have to be very careful when using :code:`FixAtPointBC` -- |
| 235 | + although similar cautionary behaviour would also have to be taken if using |
| 236 | + :code:`nullspace` and :code:`transpose_nullspace` instead. |
| 237 | + |
| 238 | +Finally, we form and solve the nonlinear variational problem for :math:`T_0 \in \{1, 10, 100, 1000, 10000 \}`:: |
| 239 | + |
| 240 | + NLVP = NonlinearVariationalProblem(F, z, bcs=[bc_u]+aux_bcs) |
| 241 | + NLVS = NonlinearVariationalSolver(NLVP) |
| 242 | + |
| 243 | + (u, p, T_aux, l) = z.subfunctions |
| 244 | + File = VTKFile(f"output/boussinesq.pvd") |
| 245 | + |
| 246 | + for i in range(0, 5): |
| 247 | + T0.assign(10**(i)) |
| 248 | + l.assign(Constant(T0)) |
| 249 | + NLVS.solve() |
| 250 | + |
| 251 | + u_out = assemble(project(u, Z.sub(0))) |
| 252 | + p_out = assemble(project(p, Z.sub(1))) |
| 253 | + T_out = assemble(project(T, Z.sub(2))) |
| 254 | + |
| 255 | + u_out.rename("u") |
| 256 | + p_out.rename("p") |
| 257 | + T_out.rename("T") |
| 258 | + |
| 259 | + File.write(u_out, p_out, T_out, time=i) |
| 260 | + |
| 261 | +The temperature and stream lines for :math:`T_0=1` and :math:`T_0=10000` are displayed below on the left and right respectively. |
| 262 | + |
| 263 | ++-------------------------+-------------------------+ |
| 264 | +| .. image:: T0_1.png | .. image:: T0_10000.png | |
| 265 | +| :width: 100% | :width: 100% | |
| 266 | ++-------------------------+-------------------------+ |
| 267 | + |
| 268 | +A Python script version of this demo can be found :demo:`here |
| 269 | +<boussinesq.py>`. |
0 commit comments