Hi,
I’m trying to run a simulation for a two-layer monoblock in 2D, but I think there are issues when applying the mesh that I have created to FESTIM when it comes to applying surface boundary conditions. When plotting the temperature field after running the simulation it does not vary as per the given boundary conditions, which are 800K on the left, top, bottom and inner curved surface, and 400K at the interface between the tungsten and CuCrZr and I think the error is with how I am defining and marking the curved surfaces.
Here is how I am generating the mesh
import numpy as np
import matplotlib.pyplot as plt
import festim as F
from fenics import RectangleMesh, MeshFunction, Point, CompiledSubDomain, plot
from scipy.constants import Avogadro
from dolfin import *
from mshr import *
left_id = 1
top_id = 2
bottom_id = 3
outer_surface_id = 4
inner_surface_id = 5
outer_layer_id = 1
inner_layer_id = 2
def create_mesh(width, height, inner_radius, outer_radius, n):
y_min = -height/2
y_max = height/2
x_max = width/2
tol = 1e-14
class InnerSemiCircleBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near((x[0]-x_max)*(x[0]-x_max)+(x[1]*x[1]), inner_radius*inner_radius, tol)
class OuterSemiCircleBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near((x[0]-x_max)*(x[0]-x_max)+(x[1]*x[1]), outer_radius*outer_radius, tol)
base = Rectangle(Point(0, y_min), Point(x_max, y_max))
inner_circle = Circle(Point(x_max, 0), inner_radius)
domain = base - inner_circle
full_mesh = generate_mesh(domain, n)
surface_markers = MeshFunction("size_t", full_mesh, full_mesh.topology().dim()-1)
surface_markers.set_all(0)
inner_circle_boundary = InnerSemiCircleBoundary()
outer_circle_boundary = OuterSemiCircleBoundary()
left_surface = CompiledSubDomain('near(x[0], 0, tol)', tol = tol)
top_surface = CompiledSubDomain('near(x[1], y_max, tol)', y_max=y_max, tol = tol)
bottom_surface = CompiledSubDomain('near(x[1], y_min, tol)', y_min = y_min, tol = tol)
left_surface.mark(surface_markers, left_id)
top_surface.mark(surface_markers, top_id)
bottom_surface.mark(surface_markers, bottom_id)
inner_circle_boundary.mark(surface_markers, inner_surface_id)
outer_circle_boundary.mark(surface_markers, outer_surface_id)
volume_markers = MeshFunction("size_t", full_mesh, full_mesh.topology().dim())
volume_markers.set_all(1)
outer_layer_volume = CompiledSubDomain('(x[0]-x_max)*(x[0]-x_max)+x[1]*x[1]>outer_radius*outer_radius+tol',
outer_radius = outer_radius,x_max = x_max, y_min = y_min, y_max = y_max, tol = tol)
inner_layer_volume = CompiledSubDomain('(x[0]-x_max)*(x[0]-x_max)+x[1]*x[1] < outer_radius*outer_radius-tol',
x_max = x_max, outer_radius=outer_radius, inner_radius = inner_radius, tol=tol)
outer_layer_volume.mark(volume_markers, outer_layer_id)
inner_layer_volume.mark(volume_markers, inner_layer_id)
return full_mesh, surface_markers, volume_markers
and then I am trying to run this through the following FESTIM simulation
def simulation(width, height, outer_radius, inner_radius, flux, mean_depth, temperatures, abs_tolerance, rel_tolerance, n, simulation_time, folder):
tungsten = F.Material(id = 1,
D_0 = 1.5e-7,
E_D = 0.265,
S_0 = 2.7e24,
E_S = 1.14,
thermal_cond = 173)
cucrzr = F.Material(id=2,
D_0=4.8e-7,
E_D = 0.42,
S_0=4.27e23,
E_S=0.39,
thermal_cond = 320)
mesh, surface_markers, volume_markers = create_mesh(width, height, inner_radius, outer_radius, n)
model = F.Simulation()
model.mesh=F.Mesh(mesh=mesh, surface_markers=surface_markers, volume_markers=volume_markers)
model.materials = [tungsten, cucrzr]
model.T = F.HeatTransferProblem(transient = False)
bcs = [F.ImplantationDirichlet(surfaces=left_id,
phi=flux,
R_p=mean_depth,
D_0=1.5e-07,
E_D=0.265),
F.DirichletBC(
surfaces = [bottom_id, top_id, inner_surface_id],
value = 0,
field = 0),
F.DirichletBC(
surfaces = [left_id, top_id, bottom_id],
value = temperatures[0],
field = 'T'),
F.DirichletBC(
surfaces=[outer_surface_id],
value = temperatures[1],
field = 'T'
),
F.DirichletBC(
surfaces = inner_surface_id,
value = temperatures[2],
field = 'T'
)]
model.boundary_conditions=bcs
model.settings = F.Settings(
absolute_tolerance=abs_tolerance,
relative_tolerance=rel_tolerance,
chemical_pot=True,
transient = False,
traps_element_type = "DG"
)
total_mobile_1 = F.TotalVolume("solute", volume=1)
total_mobile_2 = F.TotalVolume("solute", volume = 2)
retention_1 = F.TotalVolume("retention", volume =1)
retention_2 = F.TotalVolume("retention", volume = 2)
derived = [total_mobile_1, total_mobile_2, retention_1, retention_2]
model.exports = F.Exports([F.XDMFExport("T", filename = folder + "/temp.xdmf"),
F.XDMFExport("solute", filename=folder+"/solute.xdmf"),
F.XDMFExport("retention", filename=folder+"/retention.xdmf"),
F.DerivedQuantities(derived, filename = folder + "/derived_qs.csv"),
])
model.initialise()
model.run()
return model
model = simulation(15e-3, 10e-3, 4e-3, 2e-3, 1e22, 1e-8, [800,400,800], 1e12, 1e-12, 200, 1e7, "simulation")
Here is the plotting code:
from fenics import XDMFFile, FunctionSpace, Function, plot
def load_xdmf(mesh, filename, field, element="DG", counter=-1):
"""Loads a XDMF file and store its content to a fenics.Function
Args:
mesh (fenics.mesh): the mesh of the function
filename (str): the XDMF filename
field (str): the name of the field in the XDMF file
element (str, optional): Finite element of the function.
Defaults to "CG".
counter (int, optional): timestep in the file, -1 is the
last timestep. Defaults to -1.
Returns:
fenics.Function: the content of the XDMF file as a Function
"""
V = FunctionSpace(mesh, element, 1)
u = Function(V)
XDMFFile(filename).read_checkpoint(u, field, counter)
return u
mesh = model.mesh.mesh
mobile_concentration = load_xdmf(mesh, "simulation/solute.xdmf", "mobile_concentration", element = "DG")
retention = load_xdmf(mesh, "simulation/retention.xdmf", "retention")
temperature = load_xdmf(mesh, "simulation/temp.xdmf", "temperature")
CF = plot(retention)
CB = plt.colorbar(CF)
plt.show()
CF = plot(mobile_concentration)
CB = plt.colorbar(CF)
plt.show()
CF = plot(temperature)
CB = plt.colorbar(CF)
plt.show()
Any help would be appreciated
