Hi. I’m currently trying to implement the Schwarz Alternating Algorithm for an EigenValue problem. The algorithm works in a domain that has an overlap, such as:
We start by solving the eigenvalue problem on each domain separately, which has the bilinear form:
a_n(u,v) = \int_{\Omega_n} \nabla u_n \cdot \nabla v d\Omega_n -\omega^2 \int_{\Omega_n} u_n v d\Omega
With homogeneous boundary conditions in the non-overlapping boundary, and free conditions on the overlapping boundaries:
u_n = 0 \text{ on } \partial \Omega \setminus\Gamma_n
After solving for both Subdomains, we start an iterative algorithm that applies a Robin Boundary Condition in \Gamma _n , which is a function of the other subdomain. The bilinear form of this step is given by:
a_n (u_n, v) = \int_{\Omega_n} \nabla u_n \cdot \nabla v d\Omega_n -\omega^2 \int_{\Omega_n} u_n v d\Omega - \int_ {\Gamma _n} (1 - g_n) (\nabla u_n \cdot \vec{n}) v d \Gamma_n + \int_ {\Gamma_n} f_n u_n v d {\Gamma _n}
The Schwarz Alternating Algorithm is defined by:
- Do f_1 = \nabla u_2 \cdot n_2 and g_1 = u_2 |_{\Gamma _1}
- Find u_1 ^{i+1} that is solution to a_1(u_1^{i+1},v) = 0
- Do f_2 = \nabla u_1 \cdot n_1 and g_2 = u_1 |_{\Gamma _2}
- Find u_2 ^{i+1} that is solution to a_1(u_2^{i+1},v) = 0
- Until || u_1 ^{i+1} - u_2 ^{i+1} ||_{\Omega_{12}} \leq \epsilon
I’m currently having trouble in applying this Robin Boundary Condition, as I have to access information in the other Subdomain. Currently, I have this MWE to show:
from fenics import *
from dolfin import *
from petsc4py import PETSc
## Geometries and Meshes Definition
nx = 24
ny = 24
o = 0.7 # Total Overlapping Area
b = 1.5*(1-o)/2 # Total non-Overlapping Area
omega1 = RectangleMesh(Point(0.0,0.0), Point(b+o*1.5,1), nx, ny, "right/left") # First SubDomain
omega2 = RectangleMesh(Point(b,0.0), Point(1.5,1), nx, ny, "right/left") # Second SubDomain
# Function Spaces
V1 = FunctionSpace(omega1,"Lagrange",1)
V2 = FunctionSpace(omega2,"Lagrange",1)
# Defining Boundary Domains
# Omega 1
class DirichletBoundaryOmega1(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (x[0] < b+1.5*o -DOLFIN_EPS)
class RobinBoundaryOmega1(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (x[0] > b+1.5*o -DOLFIN_EPS)
# Boundary Markers
omega1_markers = MeshFunction('size_t', omega1, omega1.topology().dim()-1,0)
omega1_markers.set_all(9999)
bc_D1 = DirichletBoundaryOmega1()
bc_R1 = RobinBoundaryOmega1()
bc_D1.mark(omega1_markers, 0)
bc_R1.mark(omega1_markers, 1)
# Omega 2
class DirichletBoundaryOmega2(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (x[0] > b+DOLFIN_EPS)
class RobinBoundaryOmega2(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (x[0] < b+DOLFIN_EPS)
# Boundary Markers
omega2_markers = MeshFunction('size_t', omega2, omega2.topology().dim()-1,0)
omega2_markers.set_all(9999)
bc_D2 = DirichletBoundaryOmega2()
bc_D2.mark(omega2_markers, 0)
bc_R2 = RobinBoundaryOmega2()
bc_R2.mark(omega2_markers, 1)
## Assemblies
# Trial/Test Functions
u1 = TrialFunction(V1)
u2 = TrialFunction(V2)
v1 = TestFunction(V1)
v2 = TestFunction(V2)
# Bilinear Forms
k1 = inner(grad(u1),grad(v1))*dx
m1 = u1*v1*dx
k2 = inner(grad(u2),grad(v2))*dx
m2 = u2*v2*dx
# Matrices
# Stiffness
K1 = PETScMatrix()
K2 = PETScMatrix()
assemble(k1, tensor=K1)
assemble(k2, tensor=K2)
# Applying Boundary Conditions
bcs1 = DirichletBC(V1, 0.0, bc_D1)
bcs1.apply(K1)
bcs2 = DirichletBC(V2, 0.0, bc_D2)
bcs2.apply(K2)
# Mass
M1 = PETScMatrix()
M2 = PETScMatrix()
assemble(m1, tensor=M1)
assemble(m2, tensor=M2)
# EigenSolver
# Create the eigensolvers
omega1solver = SLEPcEigenSolver(K1,M1)
omega1solver.parameters['spectrum'] = "smallest magnitude"
omega1solver.parameters['tolerance'] = 1e-6
omega1solver.parameters['problem_type'] = "pos_gen_non_hermitian"
omega2solver = SLEPcEigenSolver(K2,M2)
omega2solver.parameters['spectrum'] = "smallest magnitude"
omega2solver.parameters['tolerance'] = 1e-6
omega2solver.parameters['problem_type'] = "pos_gen_non_hermitian"
# Solve the initial problem
omega1solver.solve(1)
omega2solver.solve(1)
# Eigenvalue Extraction
r1, c1, rx1, cx1 = omega1solver.get_eigenpair(0)
r2, c2, rx2, cx2 = omega2solver.get_eigenpair(0)
uu = Function(V)
uu.vector()[:] = rx
uu2 = Function(V2)
uu2.vector()[:] = rx2
## Schwarz Alternating Algorithm
# First Step: Do F1 = du2/dn2 G1 = u2
W2 = VectorFunctionSpace(omega2, 'P', V2.ufl_element().degree())
grad_u2 = project(grad(uu2), W2)
print(grad_u2(0.8,0.8)[0], grad_u2(0.8,0.8)[1])
print(uu2(0.8,0.8))
u1 = TrialFunction(V1)
v1 = TestFunction(V1)
dx = Measure('dx', domain = omega1)
ds = Measure('ds', subdomain_data = omega1_markers)
class OmegaF1(UserExpression):
def eval(self, value, x):
f = grad_u2(x[0], x[1])
value[0] = f[0]
value[1] = f[1]
class OmegaG1(UserExpression):
def eval(self, value, x):
g = uu2(x[0], x[1])
value[0] = g[0]
f1 = OmegaF1()
g1 = OmegaG1()
n = FacetNormal(omega1)
ak1 = inner(grad(u1),grad(v1))*dx -(1-g1)*inner(grad(u1),n)*v1*ds(1) -f1[0]*u1*v1*ds(1) -f1[1]*u1*v1*ds(1)
am1 = u1*v1*dx
K1 = PETScMatrix()
bcs1 = DirichletBC(V1, 0.0, bc_D1)
bcs1.apply(K1)
assemble(ak1, tensor=K1)
M1 = PETScMatrix()
assemble(am1, tensor=M1)
omega1solver = SLEPcEigenSolver(K1,M1)
omega1solver.parameters['spectrum'] = "smallest magnitude"
omega1solver.parameters['tolerance'] = 1e-6
omega1solver.parameters['problem_type'] = "pos_gen_non_hermitian"
omega1solver.solve(1)
r1, c1, rx1, cx1 = omega1solver.get_eigenpair(0)
print("First Iteration Omega 1 Frequency: ", sqrt(r1))
With this code, I run into the following problems:
WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
Traceback (most recent call last):
File "SCHWARZ_OverlappingMembraneEigenValueProblem.py", line 196, in <module>
ak1 = inner(grad(u1),grad(v1))*dx -(1-g1)*inner(grad(u1),n)*v1*ds(1) -f1[0]*u1*v1*ds(1) -f1[1]*u1*v1*ds(1)
File "/home/careca/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/exproperators.py", line 449, in _getitem
all_indices, slice_indices, repeated_indices = create_slice_indices(component, shape, self.ufl_free_indices)
File "/home/careca/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/index_combination_utils.py", line 163, in create_slice_indices
if int(ind) >= shape[len(all_indices)]:
IndexError: tuple index out of range
I assume it’s something with the W2 = VectorFunctionSpace function, but I can’t find the solution. Can anyone help me?