Solution discontinuous over subdomains (continuous solution needed)


53
views
0
28 days ago by

Hi everyone,

I'm working on a magnetic field simulation with the code attached below and I'm getting very strange results.

The 3D-geometry is divided into various subdomains, some of which have currents inducing the field. The solution should be a continuous function over the whole geometry, first for the vector potential A, then for its curl B. Unfortunately, that's not what I get. The PDE-solution A is computed as constant where $L\ne0$L0, and as zero everywhere else (where L = 0). So it has discontinuities on each of the subdomain boundaries. It looks as if the solution was computed for each subdomain separately, instead of connecting all the conditions to one field.

My guess would be that something is wrong with the integration measures, but I have no idea what. What's especially confusing is that the (probably) exact same approach worked fine on a 2D model. I'd really appreciate any help finding out what's wrong here. Thanks in advance.

P. S. Don't be surprised because I split up all the vectors - that's to get a better overview over the components while working and looking for the error. I've had the same problem with the vectorial formulation.

Community: FEniCS Project

Sorry, I posted before adding the code. This is it:

from __future__ import print_function
from fenics import *
from mshr import *
from math import sin, cos, pi

import matplotlib.pyplot as plt
import numpy as np

# some constants
mu_ri = 1.0e6                 
n = 20000                    

# coil currents in Ampere
I_left = 10000.0
I_right = 10000.0

# following current densities 
J_left = I_left*n/(109.4*27.88)
J_right = I_right*n/(109.4*27.88)

mesh = Mesh("Feldgenerator.xml")
subdomains = MeshFunction("size_t", mesh, "Feldgenerator_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "Feldgenerator_facet_region.xml")

V = FunctionSpace(mesh, "CG", 2)


class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        b = bool(near(x[0], -800) or near(x[0], 1200) or near(x[1], -1200) or near(x[1], 1200) or near(x[2], -800) or near(x[2], 800))
        return b

b0 = Constant(0.0)
dbc = DirichletBoundary()
bc = DirichletBC(V, b0, dbc)

dV = Measure("dx")(domain = mesh, subdomain_data = subdomains)

# define relative permeability for subdomains
class Permeability(Expression):
    def __init__(self, subdomains, **kwargs):
        self.subdomains = subdomains
    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] == 1:
            values[0] = mu_ri
        else:
            values[0] = 1.0

# initialise permeability
Mu = 4.0*pi*1.0e-7*Permeability(subdomains, degree = 0)

# set mu for whole mesh
mu = project(Mu, V)

# define current densities in x-direction for subdomains
class CurrentDensity_x(Expression):
    def __init__(self, subdomains, **kwargs):
        self.subdomains = subdomains
    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] > 1 and self.subdomains[cell.index] < 18:   
            if self.subdomains[cell.index] < 10:
                J_value = J_left
            else:
                J_value = J_right
            if (self.subdomains[cell.index] == 2 or self.subdomains[cell.index] == 14):    
                values[0] = J_value
            elif (self.subdomains[cell.index] == 3 or self.subdomains[cell.index] == 13):
                values[0] = 1/sqrt(2)*J_value
            elif (self.subdomains[cell.index] == 4 or self.subdomains[cell.index] == 12):
                values[0] = 0
            elif (self.subdomains[cell.index] == 5 or self.subdomains[cell.index] == 11):
                values[0] = -1/sqrt(2)*J_value
            elif (self.subdomains[cell.index] == 6 or self.subdomains[cell.index] == 10):
                values[0] = -J_value
            elif (self.subdomains[cell.index] == 7 or self.subdomains[cell.index] == 17):
                values[0] = -1/sqrt(2)*J_value
            elif (self.subdomains[cell.index] == 8 or self.subdomains[cell.index] == 16):
                values[0] = 0
            else:                                                                           # subdomains 9 and 15
                values[0] = 1/sqrt(2)*J_value
        else:
            values[0] = 0

# define current densities in x-direction for subdomains
class CurrentDensity_z(Expression):
    def __init__(self, subdomains, **kwargs):
        self.subdomains = subdomains
    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] > 1 and self.subdomains[cell.index] < 18:            # check if subdomain is coil part
            if self.subdomains[cell.index] < 10:
                J_value = J_left
            else:
                J_value = J_right
            if (self.subdomains[cell.index] == 2 or self.subdomains[cell.index] == 14):
                values[0] = 0
            elif (self.subdomains[cell.index] == 3 or self.subdomains[cell.index] == 13):
                values[0] = -1/sqrt(2)*J_value
            elif (self.subdomains[cell.index] == 4 or self.subdomains[cell.index] == 12):
                values[0] = -J_value
            elif (self.subdomains[cell.index] == 5 or self.subdomains[cell.index] == 11):
                values[0] = -1/sqrt(2)*J_value
            elif (self.subdomains[cell.index] == 6 or self.subdomains[cell.index] == 10):
                values[0] = 0
            elif (self.subdomains[cell.index] == 7 or self.subdomains[cell.index] == 17):
                values[0] = 1/sqrt(2)*J_value
            elif (self.subdomains[cell.index] == 8 or self.subdomains[cell.index] == 16):
                values[0] = J_value
            else:                                                                           # subdomains 9 and 15
                values[0] = 1/sqrt(2)*J_value
        else:
            values[0] = 0

# initialise current densities
JX = CurrentDensity_x(subdomains, degree = 1)
JZ = CurrentDensity_z(subdomains, degree = 1)

Jx = project(JX, V)
Jz = project(JZ, V)


# define variational problems for x, y, z
Ax = TrialFunction(V)
vx = TestFunction(V)
Ay = Function(V)
vy = TestFunction(V)
Az = TrialFunction(V)
vz = TestFunction(V)

# left hand sides of variational problem
ax = dot(grad(Ax), grad(vx))*dV
ay = dot(grad(Ay), grad(vy))*dV
az = dot(grad(Az), grad(vz))*dV

# right hand sides of variational problems
Lx = mu*Jx*vx*dV
#Ly = 0
Lz = mu*Jz*vz*dV

# solve variational problems
Ax = Function(V)
solve(ax == Lx, Ax, bc)
#Ay = Function(Vy)
solve(ay == 0, Ay, bc)
Az = Function(V)
solve(az == Lz, Az, bc)

plot(Ax)
plt.show()
plot(Ay)
plt.show()
plot(Az)
plt.show()



# compute magnetic field (B = curl(A))
W = VectorFunctionSpace(mesh, "CG", 2)
A = project(as_vector((Ax, Ay, Az)), W)

B = project(1.0e6*curl(A), W)

Bx, By, Bz = split(B)

plot(Bx)
plt.show()
plot(By)
plt.show()
plot(Bz)
plt.show()

# print solutions on screen
for v in vertices(mesh):
    x = v.point().x()
    y = v.point().y()
    z = v.point().z()
    print(x, y, z, A(Point(x, y, z))[0], A(Point(x, y, z))[1], A(Point(x, y, z))[2])
 
# save solutions A and B for plotting in paraview
vtkfile_A = File("solutions/scalar-03/potential.pvd")
vtkfile_A << A
vtkfile_B =  File("solutions/scalar-03/field.pvd")
vtkfile_B << B
written 28 days ago by Rebecca Pflanze  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »