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):
    epsilon_val = 1.0/2.0*(grad(u)+grad(u).T)
    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
com_a = rho_a/K_f*inner(af_new, Wt)*dx_a + inner(grad(w), grad(Wt))*dx_a 
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.
Thank you in advance for your help.

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 :slight_smile: Feel free to ask further questions on how to get this to work for your model