How to solve decouple problems defined in different subdomains


289
views
1
7 months ago by
Let  $\Omega=\left[0,1\right]^3$Ω=[0,1]3, $\Omega_1=\left[0.4,0.6\right]^3$Ω1=[0.4,0.6]3,  $\Omega_0=\Omega\backslash\overline{\Omega}_1$Ω0=Ω\Ω1 and  $\Gamma=\overline{\Omega}_0\cap\overline{\Omega}_1$Γ=Ω0Ω1.
Now, for  $i=1,\ldots,n$i=1,…,n, find  $\textbf{A}_i$Ai and  $\phi_i$ϕi such that
  $\int_{\Omega_1}\delta\textbf{A}_i\cdot\textbf{Q}+\int_{\Omega}\nabla\times\textbf{A}_i\cdot\nabla\times\textbf{Q}+\int_{\Omega}\nabla\cdot\textbf{A}_i\cdot\nabla\cdot\textbf{Q}=\int_{\Omega_0}\textbf{J}_i\cdot\textbf{Q}-\int_{\Omega_1}\nabla\phi_{i-1}\cdot\textbf{Q}$Ω1δAi·Q+Ω×Ai·×Q+Ω·Ai··Q=Ω0Ji·QΩ1ϕi1·Q,
  $\int_{\Omega_1}\nabla\phi_i\cdot\nabla\psi=\int_{\Omega_1}\delta\textbf{A}_i\cdot\nabla\psi$Ω1ϕi·ψ=Ω1δAi·ψ,
where
$\delta u_i=\frac{u_i-u_{i-1}}{\tau}$δui=uiui1τ .
My question is how to solve  $\phi$ϕ which is defined on  $\Omega_1$Ω1 only.
Community: FEniCS Project
I asked a similar question on the old forum. Maybe it'll help?
https://fenicsproject.org/qa/1605/2-unknowns-solved-in-2-different-subdomains
written 7 months ago by Mike Welland  
Hi Mike,
Thanks a lot for your reply. After I reading the web (and the links in it also), I tried to use DirichletBC on restrict $\phi$ϕ in $\Omega_0$Ω0, and the program can run now.
Thanks again for your help.
written 7 months ago by Osbert Wang  
As a possible solution, please see my recent post at https://www.allanswered.com/post/qvxae/announcement-multiphenics-easy-prototyping-of-multiphysics-problems-in-fenics/ . There are of course several other possible solutions, related to the one linked by Mike Welland.
written 7 months ago by Francesco Ballarin  
Hello Francesco,
I'm a amateur for FEniCS, and I think your multiphenics is a little hard to understand for me now.
Thank you all the same.
written 7 months ago by Osbert Wang  

1 Answer


1
7 months ago by

I post an example for my question here. I use DirichletBC to restrict  $\phi$ϕ .

from __future__ import print_function
from fenics import *
import numpy as np
import time 
import sys

# Begin the program
tic = time.clock()

# Define time parameters
T = 1.0
num_step = 2**3
dt = T/num_step

# Define mesh. The element type is mesh.ufl_cell() = tetrahedron
nx = ny = nz = 2**3
mesh = UnitCubeMesh(nx, ny, nz)

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 1)
W = FunctionSpace(mesh, 'P', 1)

# Define subdomain
class subdomain_1(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (0.25, 0.75)) and between(x[1], (0.25, 0.75)) and between(x[2], (0.25, 0.75)))
    
Omega_1 = subdomain_1()

# Mark subdomain
cell_marker = CellFunction('size_t', mesh)
cell_marker.set_all(0)
Omega_1.mark(cell_marker, 1)

# Redefine the integration measure on subdomain
dx = Measure('dx', domain = mesh, subdomain_data = cell_marker)

# Define interface
class interface_1(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[0], 0.25) and between(x[1], (0.25, 0.75)) and between(x[2], (0.25, 0.75)))

class interface_2(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[0], 0.75) and between(x[1], (0.25, 0.75)) and between(x[2], (0.25, 0.75)))

class interface_3(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[1], 0.25) and between(x[0], (0.25, 0.75)) and between(x[2], (0.25, 0.75)))
    
class interface_4(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[1], 0.75) and between(x[0], (0.25, 0.75)) and between(x[2], (0.25, 0.75)))
    
class interface_5(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[2], 0.25) and between(x[0], (0.25, 0.75)) and between(x[1], (0.25, 0.75)))

class interface_6(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[2], 0.75) and between(x[0], (0.25, 0.75)) and between(x[1], (0.25, 0.75)))
    
gamma_1 = interface_1()
gamma_2 = interface_2()
gamma_3 = interface_3()
gamma_4 = interface_4()
gamma_5 = interface_5()
gamma_6 = interface_6()

# Mark interface
facet_marker = FacetFunction('size_t', mesh)
facet_marker.set_all(0)
gamma_1.mark(facet_marker, 1)
gamma_2.mark(facet_marker, 2)
gamma_3.mark(facet_marker, 3)
gamma_4.mark(facet_marker, 4)
gamma_5.mark(facet_marker, 5)
gamma_6.mark(facet_marker, 6)

# Redefine the integration measure on interface. Pay attention to capital 'S'
dS = Measure('dS', domain = mesh, subdomain_data = facet_marker)

# Define exact solutions
A_e = Expression(('t*(x[0]*x[0] - 2*x[2]*x[0])', 't*(x[1]*x[1] - 2*x[0]*x[1])', 't*(x[2]*x[2] - 2*x[1]*x[2])'), t = 0.0, element = V.ufl_element())
phi_e = Expression('x[0]*x[0] + x[1]*x[1] + x[2]*x[2]', element = W.ufl_element())
E_e = Expression(('2*x[2]*x[0] - x[0]*x[0] - 2*x[0]', '2*x[0]*x[1] - x[1]*x[1] - 2*x[1]', '2*x[1]*x[2] - x[2]*x[2] - 2*x[2]'), t = 0.0, element = V.ufl_element())
B_e = Expression(('-2*t*x[2]', '-2*t*x[0]', '-2*t*x[1]'), t = 0.0, element = V.ufl_element())

# Define the Dirichlet boundary
def boundary_A(x, on_boundary):
    return on_boundary

def boundary_phi(x, on_boundary):
    return (between(x[0], (0.0, 0.25)) or between(x[0], (0.75, 1.0)) or between(x[1], (0.0, 0.25)) or between(x[1], (0.75, 1.0)) or between(x[2], (0.0, 0.25)) or between(x[2], (0.75, 1.0)))

bc_A = DirichletBC(V, A_e, boundary_A)
bc_phi = DirichletBC(W, phi_e, boundary_phi)
bc_E = DirichletBC(V, E_e, boundary_phi)

# Define material coefficient
class Conductivity(Expression):
    def __init__(self, cell_marker, sigma_0, sigma_1, **kwargs):
        self.cell_marker = cell_marker
        self.sigma_0 = sigma_0
        self.sigma_1 = sigma_1
    
    def eval_cell(self, value, x, cell):
        if self.cell_marker[cell.index] == 0:
            value[0] = self.sigma_0
        else:
            value[0] = self.sigma_1

sigma = Conductivity(cell_marker = cell_marker, sigma_0 = 0, sigma_1 = 1, element = W.ufl_element())

class Permiability(Expression):
    def __init__(self, cell_marker, mu_0, mu_1, **kwargs):
        self.cell_marker = cell_marker
        self.mu_0 = mu_0
        self.mu_1 = mu_1
        
    def eval_cell(self, value, x, cell):
        if self.cell_marker[cell.index] == 0:
            value[0] = self.mu_0
        else:
            value[0] = self.mu_1
            
mu = Permiability(cell_marker = cell_marker, mu_0 = 2, mu_1 = 1, element = W.ufl_element())

# Define current density on the whole domain
class Current_body(Expression):
    def __init__(self, cell_marker, J_0, J_1, **kwargs):
        self.cell_marker = cell_marker
        self.J_0 = J_0
        self.J_1 = J_1
        
    def eval_cell(self, values, x, cell):
        if self.cell_marker[cell.index] == 0:
            values[:] = self.J_0(x)
        else:
            values[:] = self.J_1(x)
            
    def value_shape(self): 
        return (3,) 
    
J_0 = Expression(('-t', '-t', '-t'), t = 0.0, element = V.ufl_element())
J_1 = Expression(('x[0]*x[0] - 2*x[2]*x[0] + 2*x[0] - 2*t', 'x[1]*x[1] - 2*x[0]*x[1] + 2*x[1] - 2*t', 'x[2]*x[2] - 2*x[1]*x[2] + 2*x[2] - 2*t'), t = 0.0, element = V.ufl_element())
J = Current_body(cell_marker = cell_marker, J_0 = J_0, J_1 = J_1, element = V.ufl_element())

# Define current density on the interface
class Current_surface(Expression):
    def __init__(self, facet_marker, J_s_1, J_s_2, J_s_3, J_s_4, J_s_5, J_s_6, **kwargs):
        self.facet_marker = facet_marker
        self.J_s_1 = J_s_1
        self.J_s_2 = J_s_2
        self.J_s_3 = J_s_3
        self.J_s_4 = J_s_4
        self.J_s_5 = J_s_5
        self.J_s_6 = J_s_6
        
    def eval_cell(self, values, x, facet):
        if self.facet_marker[facet.index] == 1:
             values[:] = self.J_s_1(x)
        if self.facet_marker[facet.index] == 2:
             values[:] = self.J_s_2(x)
        if self.facet_marker[facet.index] == 3:
             values[:] = self.J_s_3(x)
        if self.facet_marker[facet.index] == 4:
             values[:] = self.J_s_4(x)
        if self.facet_marker[facet.index] == 5:
             values[:] = self.J_s_5(x)
        if self.facet_marker[facet.index] == 6:
             values[:] = self.J_s_6(x)
    
    def value_shape(self):
        return (3,)

J_s_1 = Expression(('0.0*t', '-x[1]*t', 'x[0]*t'), element = V.ufl_element(), t = 0.0)
J_s_2 = Expression(('0.0*t', 'x[1]*t', '-x[0]*t'), element = V.ufl_element(), t = 0.0)
J_s_3 = Expression(('x[1]*t', '0.0*t', '-x[2]*t'), element = V.ufl_element(), t = 0.0)
J_s_4 = Expression(('-x[1]*t', '0.0*t', 'x[2]*t'), element = V.ufl_element(), t = 0.0)
J_s_5 = Expression(('-x[0]*t', 'x[2]*t', '0.0*t'), element = V.ufl_element(), t = 0.0)
J_s_6 = Expression(('x[0]*t', '-x[2]*t', '0.0*t'), element = V.ufl_element(), t = 0.0)
J_s = Current_surface(facet_marker, J_s_1, J_s_2, J_s_3, J_s_4, J_s_5, J_s_6, element = V.ufl_element())

# Define charge density on the interface
class Charge_surface(Expression):
    def __init__(self, facet_marker, rho_s_1, rho_s_2, rho_s_3, rho_s_4, rho_s_5, rho_s_6, **kwargs):
        self.facet_marker = facet_marker
        self.rho_s_1 = rho_s_1
        self.rho_s_2 = rho_s_2
        self.rho_s_3 = rho_s_3
        self.rho_s_4 = rho_s_4
        self.rho_s_5 = rho_s_5
        self.rho_s_6 = rho_s_6
        
    def eval_cell(self, value, x, facet):
        if self.facet_marker[facet.index] == 1:
             value[0] = self.rho_s_1(x)
        if self.facet_marker[facet.index] == 2:
             value[0] = self.rho_s_2(x)
        if self.facet_marker[facet.index] == 3:
             value[0] = self.rho_s_3(x)
        if self.facet_marker[facet.index] == 4:
             value[0] = self.rho_s_4(x)
        if self.facet_marker[facet.index] == 5:
             value[0] = self.rho_s_5(x)
        if self.facet_marker[facet.index] == 6:
             value[0] = self.rho_s_6(x)

rho_s_1 = Expression('2*t', t = 0.0, element = W.ufl_element())
rho_s_2 = Expression('-2*t', t = 0.0, element = W.ufl_element())
rho_s_3 = Expression('2*t', t = 0.0, element = W.ufl_element())
rho_s_4 = Expression('-2*t', t = 0.0, element = W.ufl_element())
rho_s_5 = Expression('2*t', t = 0.0, element = W.ufl_element())
rho_s_6 = Expression('-2*t', t = 0.0, element = W.ufl_element())
rho_s = Charge_surface(facet_marker, rho_s_1, rho_s_2, rho_s_3, rho_s_4, rho_s_5, rho_s_6, element = W.ufl_element())

# Define time step
tau = Constant(dt)

# Define trial and test functions
A = TrialFunction(V)
phi = TrialFunction(W)
E = TrialFunction(V)
B = TrialFunction(V)

Q = TestFunction(V)
psi = TestFunction(W)

# Define initial value
A_0 = interpolate(A_e, V)
phi_0 = interpolate(phi_e, W)

# Define the solution vector
A_1 = Function(V)
phi_1 = Function(W)
E_1 = Function(V)
B_1 = Function(V)

# Define variational formulation
a1 = dot(sigma*A, Q)*dx \
   + dot(tau/mu*curl(A), curl(Q))*dx \
   + tau/mu*div(A)*div(Q)*dx
L1 = dot(tau*J, Q)*dx \
   - dot(tau*J_s, avg(Q))*dS(1) \
   - dot(tau*J_s, avg(Q))*dS(2) \
   - dot(tau*J_s, avg(Q))*dS(3) \
   - dot(tau*J_s, avg(Q))*dS(4) \
   - dot(tau*J_s, avg(Q))*dS(5) \
   - dot(tau*J_s, avg(Q))*dS(6) \
   + dot(sigma*A_0, Q)*dx \
   - dot(tau*sigma*grad(phi_0), Q)*dx

a2 = dot(tau*sigma*grad(phi), grad(psi))*dx
L2 = dot(tau*J, grad(psi))*dx \
   - tau*rho_s*avg(psi)*dS(1) \
   - tau*rho_s*avg(psi)*dS(2) \
   - tau*rho_s*avg(psi)*dS(3) \
   - tau*rho_s*avg(psi)*dS(4) \
   - tau*rho_s*avg(psi)*dS(5) \
   - tau*rho_s*avg(psi)*dS(6) \
   - dot(sigma*A_1, grad(psi))*dx \
   + dot(sigma*A_0, grad(psi))*dx

a3 = tau*dot(E, Q)*dx
L3 = dot(A_0, Q)*dx \
   - dot(A_1, Q)*dx \
   - tau*dot(grad(phi_1), Q)*dx(1)

a4 = dot(B, Q)*dx
L4 = dot(curl(A_1), Q)*dx

# Assemble stiffness matrix
stif1 = assemble(a1)
stif2 = assemble(a2)
stif3 = assemble(a3)
stif4 = assemble(a4)

# Apply boundary conditions to stiffness matrix
bc_A.apply(stif1)
bc_phi.apply(stif2)
bc_E.apply(stif3)

# Create progress bar
progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Define vectors to record errors
error_L2_E_1 = np.zeros(num_step)
error_L2_H_1 = np.zeros(num_step)

# Time-stepping
t = 0.0
for n in range(num_step):
    # Update current time
    t += dt
    A_e.t = t
    E_e.t = t
    B_e.t = t
    J_0.t = t
    J_1.t = t
    J_s_1.t = t
    J_s_2.t = t
    J_s_3.t = t
    J_s_4.t = t
    J_s_5.t = t
    J_s_6.t = t
    rho_s_1.t = t
    rho_s_2.t = t
    rho_s_3.t = t
    rho_s_4.t = t
    rho_s_5.t = t
    rho_s_6.t = t
    
    # Assemble load vector of L1
    load1 = assemble(L1)
    
    # Apply boundary conditions to load vector of L1
    bc_A.apply(load1)
    
    # Compute solution of A
    solve(stif1, A_1.vector(), load1, 'mumps')
    
    # Assemble load vector of L2
    load2 = assemble(L2)
    
    # Apply boundary conditions to load vector of L2
    bc_phi.apply(load2)
    
    # Compute solution of phi
    solve(stif2, phi_1.vector(), load2, 'mumps')
    
    # Assemble load vector of L3
    load3 = assemble(L3)
    
    # Apply boundary conditions to load vector of L3
    bc_E.apply(load3)
    
    # Compute solution of E
    solve(stif3, E_1.vector(), load3, 'mumps')
    
    # Assemble load vector of L4
    load4 = assemble(L4)
    
    # Compute solution of B
    solve(stif4, B_1.vector(), load4, 'mumps')
    
    # Update previous solution
    phi_0.assign(phi_1)
    A_0.assign(A_1)
    
    # Update progress bar
    progress.update(t/T)
        
    # Compute L2 error on each time step
    error_L2_E_1[n] = assemble((E_e - E_1)**2*dx(1))
    error_L2_H_1[n] = assemble((B_e/mu - B_1/mu)**2*dx)

# Compute L2 error on whole time interval
error_L2_E = np.trapz(error_L2_E_1, dx = dt)
error_L2_H = np.trapz(error_L2_H_1, dx = dt)

# Compute maximum error on whole time interval
error_max_E = np.max(error_L2_E_1)
error_max_H = np.max(error_L2_H_1)

# Print errors
print('log2(error_L2_E) = %f' % np.log2(error_L2_E))
print('log2(error_L2_H) = %f' % np.log2(error_L2_H))
print('log2(error_max_E) = %f' % np.log2(error_max_E))
print('log2(error_max_H) = %f' % np.log2(error_max_H))

# End the program
toc = time.clock()
print('The totle time is %.2f' % (toc - tic))
Please login to add an answer/comment or follow this question.

Similar posts:
Search »