Using different time steps in different subdomains

variational formulation
1

Hello,

I want to implement a time-dependent problem where different time step sizes are used in subdomains. I use FEniCSx 0.4.1 but an update to the latest version would be no problem.

To illustrate what I want to achieve, consider the given picture:

image
image605Ɨ599 4.14 KB

I have a domain consisting of two subdomains. I know that due to material parameters the solution in subdomain 1 changes more within a certain period of time than in subdomain 2 (which can be considered constant within that time period). Therefore, I want to only calculate the solution in subdomain 1 for N time steps, update subdomain 2 after N time steps, and repeat. To implement this approach, I need a way to somehow fixate the solution in subdomain 2 or apply a known solution.

While searching for similar problems, I encountered this: Integrating on subdomain. But all answers there and for similar questions are given for legacy FEniCS and I struggle to apply them to FEniCSx.

As a simple model problem, I modified the heat equation form the tutorial. Now 2 different materials are used and only subdomain 1 should be calculated. The entire code is given at the end.

I came up with the following idea (assuming the subdomains are marked as 1 and 2):

a = u * v * dx(1) + dt*ufl.dot(ufl.grad(u), ufl.grad(v)) * dx(1) 
L = (u_n + dt * f) * v * dx(1)
a += u_known * dx(2)
L += u_known * dx(2)

Where I try to use a known solution in subdomain 2. Obviously the last two lines are wrong. My questions are: How can I apply a known solution in subdomain 2? Is there a different approach that I should use? Is what Iā€™m trying to do even possible in FEniCSx?

Thank you in advance!

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from dolfinx import fem, mesh, io, plot

# Define temporal parameters
t = 0 # Start time
T = 2.0 # Final time
num_steps = 61     
dt = T / num_steps # time step size

# Define mesh
nx, ny = 50, 50
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-2, -2]), np.array([2, 2])], [nx, ny], mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("CG", 1))

Q = fem.FunctionSpace(domain, ("DG", 0))

def Omega_1(x):
    return x[1] >= 0.0

def Omega_2(x):
    return x[1] <= 0.0

kappa = fem.Function(Q)
cells_1 = mesh.locate_entities(domain, domain.topology.dim, Omega_1)
cells_2 = mesh.locate_entities(domain, domain.topology.dim, Omega_2)

kappa.x.array[cells_1] = np.full_like(cells_1, 1, dtype=ScalarType)
kappa.x.array[cells_2] = np.full_like(cells_2, 0.01, dtype=ScalarType)

# Create initial condition
def initial_condition(x, a=5):
    return np.exp(-a*(x[0]**2+x[1]**2))
u_n = fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(initial_condition)

# Create boundary condition
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(
    domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
bc = fem.dirichletbc(PETSc.ScalarType(0), fem.locate_dofs_topological(V, fdim, boundary_facets), V)

xdmf = io.XDMFFile(domain.comm, "output.xdmf", "w")
xdmf.write_mesh(domain)

# Define solution variable, and interpolate initial solution for visualization in Paraview
uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)
xdmf.write_function(uh, t)

import ufl
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))

# The following two lines of code need to be changed in a way that
# only subdomain 1 is calculated and subdomain 2 remains unchanged.
a = u * v * ufl.dx + dt*kappa*ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx 
L = (u_n + dt * f) * v * ufl.dx
##################################################################

bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = fem.petsc.assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = fem.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)
    fem.petsc.assemble_vector(b, linear_form)
    
    # Apply Dirichlet boundary condition to the vector
    fem.petsc.apply_lifting(b, [bilinear_form], [[bc]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    fem.petsc.set_bc(b, [bc])

    # 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()
1 Like
2

You would need me to use a Dirichlet boundary condition on all the interior nodes of the subdomain you would like to pin for a set of time steps.

3

Thank you for your help!