Hi All,
I am working on a simulation with three different materials: Tungsten (1 cm), EURO FER (3 cm), and FLiBe (2 cm). The geometry is 2D. I successfully performed the pure diffusion problem without any issues. However, I encounter problems when I introduce velocity and the advection term in FLiBe.
I defined the velocity expression as follows:
from fenics import interpolate, Expression, VectorFunctionSpace
functionspace = VectorFunctionSpace(my_model_adv.mesh.mesh, "CG", 1)
velocity = Expression(
("0", "(x[0]>=0.4 && x[0]<=0.6) ? 2:0"),
degree=0
)
velocity = interpolate(velocity, functionspace)
CS = plot(velocity, title="velocity")
plt.colorbar(CS)
plt.show()
However, the plot is completely empty and the colorbar is meaningless. If I substitute the expression from the workshop Task 06 , the vectors are displayed correctly.
I also have another problem when trying to export the concentration using:
my_model_adv.exports = [F.XDMFExport("solute", folder=results_folder)]
It seems that the solution is not updated, and it exports the solution from the pure diffusion case. I do not understand why this happens.
Below, you can find the full code for reference:
import numpy as np
import meshio
import matplotlib.pyplot as plt
import festim as F
import os
import copy
from fenics import plot
msh = meshio.read("/home/ndrea_i_aio/festim/TESI/Mesh_Tun_EURO_FLiBe_NEW.med")
print("CELLS:", msh.cells)
print("CELL DATA:", msh.cell_data_dict.keys())
print("CELL SETS:", msh.cell_sets)
print("FIELD DATA:", msh.field_data)
def convert_med_to_xdmf(
med_file,
cell_file="mesh_domains.xdmf",
facet_file="mesh_boundaries.xdmf",
cell_type="triangle",
facet_type="line",
):
"""Converts a MED mesh to XDMF
Args:
med_file (str): the name of the MED file
cell_file (str, optional): the name of the file containing the
volume markers. Defaults to "mesh_domains.xdmf".
facet_file (str, optional): the name of the file containing the
surface markers.. Defaults to "mesh_boundaries.xdmf".
cell_type (str, optional): The topology of the cells. Defaults to "tetra".
facet_type (str, optional): The topology of the facets. Defaults to "triangle".
Returns:
dict, dict: the correspondance dict, the cell types
"""
msh = meshio.read(med_file)
correspondance_dict = msh.cell_tags
cell_data_types = msh.cell_data_dict["cell_tags"].keys()
for mesh_block in msh.cells:
if mesh_block.type == cell_type:
meshio.write_points_cells(
cell_file,
msh.points,
[mesh_block],
cell_data={"f": [-1 * msh.cell_data_dict["cell_tags"][cell_type]]},
)
elif mesh_block.type == facet_type:
meshio.write_points_cells(
facet_file,
msh.points,
[mesh_block],
cell_data={"f": [-1 * msh.cell_data_dict["cell_tags"][facet_type]]},
)
return correspondance_dict, cell_data_types
correspondance_dict, cell_data_types = convert_med_to_xdmf("TESI/Mesh_Tun_EURO_FLiBe_NEW.med", cell_file="TESI/mesh_domains.xdmf", facet_file="TESI/mesh_boundaries.xdmf")
mesh = meshio.read("TESI/mesh_domains.xdmf")
print(mesh)
mesh_bnd = meshio.read("TESI/mesh_boundaries.xdmf")
print(mesh_bnd)
"""
# Carica dominio (triangoli)
domain = meshio.read("TESI/mesh_domains.xdmf")
# Carica bordi (linee)
boundary = meshio.read("TESI/mesh_boundaries.xdmf")
# Punti della mesh
points = domain.points
# --- Plotta triangoli (dominio)
triangles = domain.cells_dict.get("triangle")
if triangles is not None:
plt.triplot(points[:, 0], points[:, 1], triangles, color="lightblue")
# --- Plotta linee (bordi)
lines = boundary.cells_dict.get("line")
if lines is not None:
for line in lines:
x = points[line, 0]
y = points[line, 1]
plt.plot(x, y, color="red", linewidth=1.5)
plt.gca().set_aspect("equal")
plt.title("Mesh check: dominio (blu) + bordi (rosso)")
plt.show()
print(correspondance_dict)
"""
my_model = F.Simulation()
my_model.mesh = F.MeshFromXDMF(volume_file="TESI/mesh_domains.xdmf", boundary_file="TESI/mesh_boundaries.xdmf")
tungsten = F.Material(
id=12,
D_0=4.1e-7/(np.sqrt(3)),
E_D=0.39,
S_0=1.35e22,
E_S=0.29,
)
euro_fer = F.Material(
id=13,
D_0=4.57e-7/(np.sqrt(3)),
E_D=0.23,
S_0=1.35e22,
E_S=0.16,
)
FLiBe = F.Material(
id=14,
D_0=9.3e-7,
E_D=42.5/96.48,
S_0=(7.9e-2)*(6.022e23),
E_S=35/96.48,
)
my_model.materials = [tungsten, euro_fer, FLiBe]
my_model.T= 676
implatantion_flux=F.ImplantationDirichlet(surfaces=7,phi=7.15e12, R_p=9.52e-10, D_0=tungsten.D_0, E_D=tungsten.E_D)
# ho un dubbio du phi dovrebbe essere atomi/m^2/s invece io ho atomi/m^2
# va bene la distanza di impiantazione?
Neumann_top_bottom_outlet=F.FluxBC(surfaces=[11, 6, 8], value=0, field=0)
Diriclet_con=F.DirichletBC(surfaces=10, value=4.6e-5*(6.022e23), field=0)
my_model.boundary_conditions=[implatantion_flux, Neumann_top_bottom_outlet, Diriclet_con]
my_model.sources=[F.Source(volume=14, field=0, value=1.1e19)]
my_model.settings=F.Settings(
absolute_tolerance=1e-11,
relative_tolerance=1e-11,
transient=False,
#final_time=1
)
#my_model.dt= F.Stepsize(initial_value=1/20)
results_folder = "TESI"
my_model.exports = [F.XDMFExport("solute", folder=results_folder)]
my_model.initialise()
my_model.run()
hydrogen_concentration = my_model.h_transport_problem.mobile.solution
cs=plot(hydrogen_concentration, title=" H transport mobile pure diffusion")
plt.colorbar(cs)
plt.show()
#plt.savefig(os.path.join(results_folder, "h_transport_problem.mobile.solution.png")) # salva immagine su file
##### ADVECTION #########
my_model_adv = F.Simulation()
my_model_adv.mesh = my_model.mesh
my_model_adv.materials = my_model.materials
my_model_adv.boundary_conditions = my_model.boundary_conditions
my_model_adv.sources = my_model.sources
my_model_adv.settings = my_model.settings
my_model_adv.T= 676
#my_model_adv.dt=my_model.dt
from fenics import interpolate, Expression, VectorFunctionSpace
functionspace = VectorFunctionSpace(my_model_adv.mesh.mesh, "CG", 1)
velocity = Expression(
("0", "(x[0]>=0.4 && x[0]<=0.6) ? 2:0"),
degree=0
)
velocity = interpolate(velocity, functionspace)
CS = plot(velocity, title="velocity")
plt.colorbar(CS)
plt.show()
from fenics import inner, dot, grad
my_model_adv.initialise()
results_folder = "TESI/adv"
my_model_adv.exports = [F.XDMFExport("solute", folder=results_folder)]
hydrogen_concentration = my_model_adv.h_transport_problem.mobile.solution
test_function_solute = my_model_adv.h_transport_problem.mobile.test_function
advection_term = inner(dot(grad(hydrogen_concentration), velocity), test_function_solute) * my_model.mesh.dx(14)
my_model_adv.h_transport_problem.F += advection_term
print(my_model_adv.h_transport_problem.F)
my_model_adv.run()
cp=plot(hydrogen_concentration,title=" H transport mobile advection diffusion")
plt.colorbar(cp)
#plot(velocity)
plt.show()
Could you please advise on how to resolve these issues?
Thank you in advance for your help.
Best regards,
