Heat equation in 3D with discontinuous source term

63
views
0
8 weeks ago by
I am trying to solve the heat equation for a vector valued quantity $\vec A (\vec r, t)$. The problem with open boundary conditions is given as
$\frac{\partial \vec A}{\partial t} - \Delta \vec A = \vec j \\ \vec A =0 ~~~~~~ \text{on \partial \Omega \times [0,T]} \\ \nabla A_i =0 ~~~~~~ \text{on \partial \Omega \times [0,T] for i =1,2,3}$
The domain of computation $\Omega$ is build from an "inner" mesh $\Omega_1$ and an "external" cubic mesh around it such that $\Omega$=$\Omega_1 \cup \Omega_2$. The source term $\vec j$ represents a current distribution and is defined as

$\vec j=\begin{cases} \vec j_0(\vec r, t), & \text{in} ~ \Omega_1\\ 0, & \text{in} ~ \Omega_2 \end{cases}$

with $\vec j_0 \cdot \vec n = 0$ on $\partial \Omega_1$. Currently I am using first order Lagrange elements for $\vec A$ and the solution seems to be ok. However, it is known that the parallel component of the magnetic field $\vec B \times \vec n$=$(\nabla \times \vec A)\times \vec n$ is in general a discontinuous function across the interface of two domains if surface currents are present. Therefore I wonder if Lagrange elements are a good choice for $\vec A$ and if it is not better to use edge elements.
Another question is if the efficiency of my code can be improved by performing the calculations in parallel. It seems to be an interesting idea to split the problem into three independent equations for the components  $A_i$ of $\vec A$ and solve the resulting system either with OpenMPI or the module multiprocessing. Unfortunately I did not find an easily understandable tutorial how to solve several equations simultaneously with fenics. I hope someone can help me with this.

File attached: mesh3.xml (813.73 KB)

File attached: mesh3_physical_region.xml (312.7 KB)

from dolfin import *
import numpy as np
import scipy
import math

import sympy as sp
from sympy.abc import a, b, c

#Constants==========================================================================================================================

dt = 1E-2
T  = 1E-2

#mesh===============================================================================================================================

mesh = Mesh('mesh3.xml')
domain_markers = MeshFunction('size_t', mesh, 'mesh3_physical_region.xml')

torus = SubMesh(mesh, domain_markers, 0)

#Function spaces and trial/test functions=============================================================================================

Ws       =       FunctionSpace(mesh, 'CG', 1)
Wv       = VectorFunctionSpace(mesh, 'CG', 1)
V_tor    =       FunctionSpace(torus,'CG', 1)

V1     = FunctionSpace(mesh, 'CG', 1)
a1_tr  = TrialFunction(V1)
a1_te  = TestFunction(V1)
a1     = Function(V1)

V2     = FunctionSpace(mesh, 'CG', 1)
a2_tr  = TrialFunction(V2)
a2_te  = TestFunction(V2)
a2     = Function(V2)

V3     = FunctionSpace(mesh, 'CG', 1)
a3_tr  = TrialFunction(V3)
a3_te  = TestFunction(V3)
a3     = Function(V3)

dxc   = Measure('dx', domain=mesh, subdomain_data=domain_markers)

#functions===============================================================================================================================

def boundary(x,on_boundary):
return on_boundary

def make_submesh_local_to_global_map(V, Vr):
u                   = Function(V)
u.set_allow_extrapolation(True)
u.vector()[:]       = np.array(range(V.dim()))
ur                  = interpolate(u, Vr)
ur.set_allow_extrapolation(True)
local_to_global_map = np.round(ur.vector().array()).astype(int)
return local_to_global_map

#Current==============================================================================================================================

x, y= sp.symbols('x[0], x[1]')

j_x1_ = -y
j_y1_ =  x
j_pre1 = 1/sp.sqrt(x**2+y**2+0.0000000001)*600

j_x1 = j_x1_*j_pre1
j_x_code = sp.printing.ccode(j_x1)
jx_exp = Expression(j_x_code, degree= 3)

j_y1 = j_y1_*j_pre1
j_y_code = sp.printing.ccode(j_y1)
jy_exp = Expression(j_y_code, degree= 3)

jx_u_ = interpolate(jx_exp, V_tor)
jy_u_ = interpolate(jy_exp, V_tor)

torus_mesh  =  make_submesh_local_to_global_map(Ws, V_tor)

jx_u1 = Function(Ws)
jx_u1.vector()[torus_mesh] = jx_u_.vector().array()
jx_u=interpolate(jx_u1, Ws)

jy_u1 = Function(Ws)
jy_u1.vector()[torus_mesh] = jy_u_.vector().array()
jy_u=interpolate(jy_u1, Ws)

#initial functions and boundary conditions=================================================================================================

a1h0 = Function(V1)
a2h0 = Function(V2)
a3h0 = Function(V3)

zeros_1 = Constant(0)
bc1     = DirichletBC(V1, zeros_1, boundary)
bc2     = DirichletBC(V2, zeros_1, boundary)
bc3     = DirichletBC(V3, zeros_1, boundary)

#Solver===============================================================================================================================

solverA = KrylovSolver("bicgstab",'ilu')
solverA.parameters["absolute_tolerance"] = 1E-15
solverA.parameters["relative_tolerance"] = 1E-15
solverA.parameters["maximum_iterations"] = 10000
solverA.parameters["nonzero_initial_guess"] = True
solverA.parameters['monitor_convergence'] = True

#bilinear forms===============================================================================================================================

AA1 = assemble(aa1)

AA2 = assemble(aa2)

AA3 = assemble(aa3)

#linear forms====================================================================================================================================

def solver_A(a, u, bc, result, kind):
if kind == 1:
b = jx_u*a1_te*dxc(0) + (1/dt)*a1h0*a1_te*(dxc(0)+dxc(1))
if kind == 2:
b = jy_u*a1_te*dxc(0) + (1/dt)*a2h0*a1_te*(dxc(0)+dxc(1))
if kind == 3:
b = (1/dt)*a3h0*a1_te*(dxc(0)+dxc(1))
b1 = assemble(b)
bc.apply(a, b1)
solverA.solve(a, u.vector(), b1)
assign(result, u)
return result

nn = 1
tc = dt
while (tc < T + 1e-10):

#------------------------------------------------------------------------------------------
print("solve A1")
a1h0_ = solver_A(AA1,a1,bc1, a1h0, 1)
#------------------------------------------------------------------------------------------

#------------------------------------------------------------------------------------------
print("solve A2")
a2h0_ = solver_A(AA2,a2,bc2, a2h0, 2)
#------------------------------------------------------------------------------------------

#------------------------------------------------------------------------------------------
print("solve A3")
a3h0_ = solver_A(AA3,a3,bc3, a3h0, 3)
#------------------------------------------------------------------------------------------

print( "TC = %12.5f" % tc)

tc = tc + dt

nn=nn+1​
Community: FEniCS Project