6.3. Time derivatives of fields on moving meshes
We will now move a mesh around and solve a diffusion equation on the moving mesh to see what happens. First of all, the mesh shall now move to the left and right by a sinus oscillation, which can be realized as follows
from pyoomph import *
from pyoomph.expressions import *
class MovingMesh(Equations):
def define_fields(self):
self.activate_coordinates_as_dofs()
def define_residuals(self):
x,xtest=var_and_test("mesh_x")
xi,t=var(["lagrangian_x","time"])
desired_pos=xi+0.125*sin(2*pi*t)
self.add_residual(weak(x-desired_pos,xtest,lagrangian=True) )
We have omitted the coordinate_space argument in activate_coordinates_as_dofs(), so that the coordinate space is automatically adjusted by the space of the diffusion equation. The diffusion equation with diffusivity \(0.01\) could be written as follows:
class DiffusionEquation(Equations):
def define_fields(self):
self.define_scalar_field("c","C2")
def define_residuals(self):
c,ctest=var_and_test("c")
# Note the ALE=False statement here. We come to it in the text
self.add_residual(weak(partial_t(c,ALE=False),ctest)+weak(0.01*grad(c),grad(ctest)))
And the Problem class combining the equations reads
class ALEProblem(Problem):
def define_problem(self):
self.add_mesh(RectangularQuadMesh(N=32,lower_left=[-0.5,-0.5]))
eqs=MovingMesh()
eqs+=MeshFileOutput()
eqs+=DiffusionEquation()
eqs+=DirichletBC(mesh_y=True)
x=var("coordinate")
eqs+=InitialCondition(c=exp(-dot(x,x)*100))
self.add_equations(eqs@"domain")
if __name__=="__main__":
with ALEProblem() as problem:
problem.run(1,numouts=20)
Note how we pin the \(y\)-coordinate by a DirichletBC applied on the entire "domain". There is no equation for the \(y\)-coordinate, so we have to fix it. We furthermore set a Gaussian initial condition for the diffusion field.
Fig. 6.2 When the mesh is moving oscillatory to the left and right, all fields move along with it. The center of the Gaussian spot is always in the center of the mesh, not of the Eulerian coordinate system. When the keyword argument ALE of the function partial_t() is set to True or "auto", it is compensated for the mesh motion. Thereby, the maximum of the field stays in the center of the coordinates, not of the mesh. Since it gets into contact with the boundaries, the field is slightly deformed.
What happens is naively counter-intuitive, namely the diffusion field \(c\) is moving along with the oscillating mesh, see Fig. 6.2. One might expect, that the maximum of the field \(c\) will be always at \(\vec{x}=0\), but it is not the case. Instead, it will be always in the center of the mesh, i.e. at \(\vec{\xi}=0\). This can be understood, when considering that fields are always approximated as functions of the nodal values \(c_l(t)\). Here, we have
where \(\psi(\vec{x},t)\) are the shape/basis functions, which are due to the moving mesh now also a function of time \(t\) and \(l\) is a summation over all nodes. When adding ALE=False to partial_t(c), we will just calculate the temporal derivatives of the coefficients \(c_l(t)\), i.e.
Thereby, when the mesh moves, the entire field (including the time derivative) will co-move with the mesh. If we want to compensate for the mesh motion, we have to compensate for the term originating from the chain rule due to the time-dependence of the mesh coordinates. To that end partial_t() has an optional argument ALE, which defaults to ALE="auto". If ALE is False, we calculate the time derivative according to (6.2). However, if ALE=True is passed, we get
Thereby, the field \(c\) is moved against the mesh motion and hence stays in place when the mesh moves. Since ALE=True requires the evaluation of the mesh velocity \(\dot{\vec{x}}\), it is not required on static meshes. On a static mesh, \(\dot{\vec{x}}=0\) holds, and hence the calculation is redundantly time-consuming during the assembly of the system. Pyoomph also allows to pass ALE="auto" (the default value) to set it to False if the mesh is static, i.e. no Equations are added that call activate_coordinates_as_dofs(). Thereby, the redundant computation of \(\dot{\vec{x}}\) is not carried out. If the mesh is moving, i.e. an equation for the mesh coordinates is present, ALE="auto" will become ALE=True, i.e. expanding according to (6.3).
Warning
The predefined variables var("coordinate") and var("mesh") are in principle the same (and so there components var("coordinate_x"), var("mesh_x"), etc). However, there are two fundamental differences: var("mesh") can have test functions, when activate_coordinates_as_dofs() has been called. Furthermore, partial_t(var("mesh")) may be non-zero on a moving mesh, i.e. yielding the mesh velocity \(\dot{\vec{x}}\), whereas partial_t(var("coordinate")) is always zero.
In conclusion, if one has a moving mesh, but want to keep spatio-temporal fields independent of the mesh motion, one should augment all partial_t() calls with an ALE="auto" (or leave it out, since it is the default value). This has been done in the bottom row of Fig. 6.2, where the weak formulation of the diffusion equation has been changed to
self.add_residual(weak(partial_t(c,ALE="auto"),ctest)+weak(0.01*grad(c),grad(ctest)))
The field \(c\) is now not following the oscillatory motion of the mesh, but stays in place.