Hello,
I’m calculating the heat-problem from the tutorial. I am using an explicit Euler method in the bottom half of the domain and an implicit Euler method in the top half of the domain.
# Variational formulation using an explicit time-discretisation
a = u * v * dx(1)
L = (u_n + dt * f) * v * dx(1) - dt * dot(grad(u_n), grad(v)) * dx(1)
L += -dt * (dot(jump(u_n, n), avg(grad(v))) - dot(avg(grad(u_n)), jump(v, n))) * dS
L += -dt * (dot(u_n * n, grad(v)) - dot(grad(u_n), v * n)) * ds
L += -dt * dot(jump(u_n, n), jump(v, n)) * dS
L += -dt * dot(u_n * n, v*n) * ds
# Variational formulation using an implicit time-discretisation
a += u * v * dx(2) + dt * dot(grad(u), grad(v)) * dx(2)
a += dt * (dot(jump(u, n), avg(grad(v))) - dot(avg(grad(u)), jump(v, n))) * dS
a += dt * (dot(u * n, grad(v)) - dot(grad(u), v * n)) * ds
a += dt * dot(jump(u, n), jump(v, n)) * dS
a += dt * dot(u * n, v * n) * ds
L += (u_n + dt * f) * v * dx(2)
The integration-measures dx(1) and dx(2) are restricted to their respective subdomains. But the boundary terms occur twice, in the linear form and in the bilinear form. The result seems to be correct. If the number of time steps is reduced enough, the subdomain with the explicit Euler method becomes unstable first which indicates that FEniCS is indeed using two different methods.
My question is: What exactly is happening in the background in this example. I can not really explain what FEniCS is doing if I provide a variational formulation like that.
Thank you!
Here is the whole code:
import gmsh
import numpy as np
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar, form, locate_dofs_geometrical, locate_dofs_topological, petsc)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities
from dolfinx.plot import create_vtk_mesh
from ufl import (SpatialCoordinate, TestFunction, TrialFunction, FacetNormal, dx, grad, inner, Measure, dot, jump, avg)
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
# Define temporal parameters
t = 0 # Start time
T = 0.1 # Final time
num_steps = 500
dt = T / num_steps # time step size
# The following code is copied from the FEniCS Tutorial #################################
gmsh.initialize()
proc = MPI.COMM_WORLD.rank
top_marker = 2
bottom_marker = 1
left_marker = 1
if proc == 0:
# We create one rectangle for each subdomain
gmsh.model.occ.addRectangle(0, 0, 0, 1, 0.5, tag=1)
gmsh.model.occ.addRectangle(0, 0.5, 0, 1, 0.5, tag=2)
# We fuse the two rectangles and keep the interface between them
gmsh.model.occ.fragment([(2,1)],[(2,2)])
gmsh.model.occ.synchronize()
# Mark the top (2) and bottom (1) rectangle
top, bottom = None, None
for surface in gmsh.model.getEntities(dim=2):
com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
if np.allclose(com, [0.5,0.25, 0]):
bottom = surface[1]
else:
top = surface[1]
gmsh.model.addPhysicalGroup(2, [bottom], bottom_marker)
gmsh.model.addPhysicalGroup(2, [top], top_marker)
# Tag the left boundary
left = []
for line in gmsh.model.getEntities(dim=1):
com = gmsh.model.occ.getCenterOfMass(line[0], line[1])
if np.isclose(com[0], 0):
left.append(line[1])
gmsh.model.addPhysicalGroup(1, left, left_marker)
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
gmsh.finalize()
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh("mesh.msh", MPI.COMM_WORLD, gdim=2)
import meshio
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
if proc == 0:
# Read in mesh
msh = meshio.read("mesh.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mt.xdmf", line_mesh)
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
ft = xdmf.read_meshtags(mesh, name="Grid")
#########################################################################################
V = FunctionSpace(domain, ("DG", 1))
# Create initial condition
def initial_condition(x, a=25):
return np.exp(-a*((x[0] - 0.5)**2+(x[1] - 0.5)**2))
u_n = Function(V)
u_n.name = "u_n"
u_n.interpolate(initial_condition)
from dolfinx import io
with io.XDMFFile(domain.comm, "output.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
dx = Measure("dx", domain, subdomain_data=ct)
ds = Measure("ds", domain, subdomain_data=ft)
dS = Measure("dS", domain, subdomain_data=ft)
# Define solution variable, and interpolate initial solution for visualization in Paraview
uh = Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)
xdmf.write_function(uh, t)
u, v = TrialFunction(V), TestFunction(V)
f = Constant(domain, ScalarType(0))
n = FacetNormal(domain)
# Variational formulation using an explicit time-discretisation
a = u * v * dx(1)
L = (u_n + dt * f) * v * dx(1) - dt * dot(grad(u_n), grad(v)) * dx(1)
L += -dt * (dot(jump(u_n, n), avg(grad(v))) - dot(avg(grad(u_n)), jump(v, n))) * dS
L += -dt * (dot(u_n * n, grad(v)) - dot(grad(u_n), v * n)) * ds
L += -dt * dot(jump(u_n, n), jump(v, n)) * dS
L += -dt * dot(u_n * n, v*n) * ds
# Variational formulation using an implicit time-discretisation
a += u * v * dx(2) + dt * dot(grad(u), grad(v)) * dx(2)
a += dt * (dot(jump(u, n), avg(grad(v))) - dot(avg(grad(u)), jump(v, n))) * dS
a += dt * (dot(u * n, grad(v)) - dot(grad(u), v * n)) * ds
a += dt * dot(jump(u, n), jump(v, n)) * dS
a += dt * dot(u * n, v * n) * ds
L += (u_n + dt * f) * v * dx(2)
bilinear_form = form(a)
linear_form = form(L)
A = petsc.assemble_matrix(bilinear_form)
A.assemble()
b = petsc.create_vector(linear_form)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
for i in range(num_steps):
t += dt
# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
loc_b.set(0)
petsc.assemble_vector(b, linear_form)
# Solve linear problem
solver.solve(b, uh.vector)
uh.x.scatter_forward()
# Update solution at previous time step (u_n)
u_n.x.array[:] = uh.x.array
# Write solution to file
xdmf.write_function(uh, t)
xdmf.close()