How to handle the acoustic-elastic interface in dolfin

Dear all,

Thinking about a problem of wave propagation through the acoustic-elastic interface I’m trying to solve the acoustic equation and elastic equation on two coupled domains, respectively. The problem is described as follows:

\rho_e \ddot{u} = \nabla \cdot \sigma on \Omega_e
\rho_f \ddot{w} = K_f \nabla^2 w on \Omega_a

subject to interface and coupling conditions:

\sigma_e\cdot n_e = -\nabla w \cdot n_a

where the n_e and n_a are outward normal vectors on the boundaries or the interface.

I defined the two subdomains and manage a MixfunctionSpace to solve it. however, I can’t get the linear system using lhs and rhs. And I don’t how to define the interface condition.

Here is my code:

from dolfin import *
import numpy as np
from scipy import signal
from mpi4py import MPI
comm = MPI.COMM_WORLD
#geometry
x0 =0.0
x1 = 1000
y0 = 0.0
y1 = 800
nx = 250
ny = 200
# Define mesh
mesh = RectangleMesh(comm, Point(x0, y0), Point(x1, y1), nx, ny, 'crossed')

## define the submeshs
domain = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
for c in cells(mesh):
domain[c] = c.midpoint().y() < 400

marker = MeshFunction("size_t", mesh, mesh.topology().dim()-1 , 0)
for f in facets(mesh):
marker[f] = 400 - 1- DOLFIN_EPS < f.midpoint().y() < 400 + 1 + DOLFIN_EPS

# Meshes
mesh_a = MeshView.create(domain, 0)
mesh_e = MeshView.create(domain, 1)
mesh_I = MeshView.create(marker, 1)

#File("domain.pvd") << domain
#File("facts.pvd") << marker
##define the functionspace
V_a = FunctionSpace(mesh_a, "CG", 1)
V_e = VectorFunctionSpace(mesh_e, "CG", 1)
V_I = FunctionSpace (mesh_I, "CG", 1)

W = MixedFunctionSpace(V_a, V_e, V_I)
(w, u, p) = TrialFunctions(W)
(Wt, Ut, Pt) = TestFunctions(W)

##Subdomain1
w_0 = Constant(0)
w_n = interpolate(w_0,V_a)
vf_n = Function(V_a)
af_n = Function(V_a)

#Subdomain2
u_0 = Constant((0.,0.))
u_n = interpolate(u_0,V_e)
vs_n = Function(V_e)
as_n = Function(V_e)

# Measures
dx_a = Measure("dx", domain=W.sub_space(0).mesh())
ds_a = Measure("ds", domain=W.sub_space(0).mesh())
dx_e = Measure("dx", domain=W.sub_space(1).mesh())
ds_e = Measure("ds", domain=W.sub_space(1).mesh())
dx_I = Measure("dx", domain=W.sub_space(2).mesh())

##parameters
rho_a = 1000
K_f = 2.25e9
rho_e = 2650
mu_e = 3.0e9
lambda_e = 4.7e9

# Stress
def sigma_e(u):
sigma_e_val = 2.0*mu_e*epsilon(u) + lambda_e*div(u)*Identity(len(u))
return sigma_e_val

# Strain
def epsilon(u):
return epsilon_val

# Time-stepping parameters
T       = 0.6
nsteps  = 1200
dt = T/nsteps
#n1 = FacetNormal(submesh1)
#n2 = FacetNormal(submesh2)
gamma   = Constant(0.5)
beta    = Constant((gamma+0.5)**2/4.)
theta   = Constant(0.5)

def update_a(u, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
beta_ = beta
else:
dt_ = float(dt)
beta_ = float(beta)
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old

def update_v(a, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
gamma_ = gamma
else:
dt_ = float(dt)
gamma_ = float(gamma)
return v_old + dt_*((1-gamma_)*a_old + gamma_*a)

def update_fields(u, u_n, vs_n, as_n):

# Update solution at last time step
# Get vectors (references)
u_vec, u0_vec  = u.vector(), u_n.vector()
vs0_vec, as0_vec = vs_n.vector(), as_n.vector()

# use update functions using vector arguments
as_vec = update_a(u_vec, u0_vec, vs0_vec, as0_vec)
vs_vec = update_v(as_vec, u0_vec, vs0_vec, as0_vec)

# Update (u_old <- u)
vs_n.vector()[:], as_n.vector()[:] = vs_vec, as_vec
u_n.vector()[:] = u.vector()

af_new = update_a(w, w_n, vf_n, af_n, ufl=True)
vf_new = update_v(af_new, w_n, vf_n, af_n, ufl=True)

as_new = update_a(u, u_n, vs_n, as_n, ufl=True)
vs_new = update_v(as_new, u_n, vs_n, as_n, ufl=True)

bcs = []
##Subdomain1
a_a, L_a = lhs(com_a), rhs(com_a)

##Subdomain2
com_e = rho_e*inner(as_new, Ut)*dx_e + inner(sigma_e(u), sym(grad(Ut)))*dx_e
a_e, L_e = lhs(com_e), rhs(com_e)



here I get the error message compute_form_with_arity cannot handle parts.

about the interface coupling, I want to write it as

n_a = FaceNormal(mesh_a)
n_e = FaceNormal(mesh_e)
a = a_a + a_e \
+ inner(p,  gard(Wt)*n_a + sigma_e(Ut)*n_e)*dI + inner(sigma_e(u)*n_e + grad(w)*n_a, Pt)*dI
L = L_a + L_e


sorry for too many lines as it is time-dependent problem, I used the Newmark method.

Hi,

As far as I know, fenics-mixed-dim does not support the use of functions like assemble(a), assemble(L), rhs and lhs. A bit of care is needed here; sometimes these functions run without raising errors but return the wrong system.

The safest option is to assemble and solve the entire system inside the time loop, using

solve(a == L, sol, solver_parameters={"linear_solver":"direct"})


Depending on your model, this might be prohibitively slow due to the reassembly of the full system.

To assemble a and L outside of the time loop, you can do a bit of manual wrangling

system = assemble_mixed_system(a == L, sol)
matrix_blocks = system[0]
rhs_blocks = system[1]

alist = extract_blocks(a)
Llist = extract_blocks(L)
matrix_blocks = [assemble_mixed(a) for a in alist]

for i in range(0, n):
rhs_blocks = [assemble_mixed(L) for L in Llist ]
matrix_blocks[2] = assemble_mixed(alist[2]) # reassemble time-dep blocks in a
matrix_blocks[5] = assemble_mixed(alist[5])

# solve
A = PETScNestMatrix(matrix_blocks)
b = Vector()
A.init_vectors(b, rhs_blocks)
A.convert_to_aij()

x = Vector()
solve(A,x,b)


I haven’t tested this on your specific example, but hopefully this helps you get started Feel free to ask further questions on how to get this to work for your model