Hi all,
I’m trying to run a steady-state, advection only simulation using the following code. I have set FixedConcentrationBC of 1e15 at the inlet, but after running the simualtion. c_in = F.AverageSurface(species, inlet) is returning as 0.
Any help would be much appreciated.
Below is my code:
import numpy as np
from mpi4py import MPI
from dolfinx.mesh import create_mesh, meshtags, exterior_facet_indices
from scipy.spatial import cKDTree
import pyvista
import basix
import ufl
from dolfinx.fem import functionspace, Function
from dolfinx.fem import assemble_scalar, form
from dolfinx.fem import locate_dofs_topological
from foam2dolfinx import OpenFOAMReader
def tag_boundary_patch(dolfinx_mesh, patch_dataset, patch_id, tol=1e-6):
fdim = dolfinx_mesh.topology.dim - 1
dolfinx_mesh.topology.create_connectivity(fdim, 0)
dolfinx_mesh.topology.create_connectivity(0, fdim)
dolfinx_mesh.topology.create_connectivity(fdim, dolfinx_mesh.topology.dim)
facet_indices = exterior_facet_indices(dolfinx_mesh.topology)
x = dolfinx_mesh.geometry.x
patch_points = patch_dataset.points
tree = cKDTree(x)
matched_vertex_indices = tree.query_ball_point(patch_points, tol)
matched_vertex_indices = list(set(i for sub in matched_vertex_indices for i in sub))
matched_facets = []
for facet in facet_indices:
vertices = dolfinx_mesh.topology.connectivity(fdim, 0).links(facet)
if all(v in matched_vertex_indices for v in vertices):
matched_facets.append(facet)
print(f"Tagging {len(matched_facets)} facets for patch ID {patch_id}")
return np.array(matched_facets, dtype=np.int32), np.full(len(matched_facets), patch_id, dtype=np.int32)
import festim as F
import h_transport_materials as htm
import ufl
from mpi4py import MPI
def FESTIM(T, time):
# Read Mesh and Velocity Fields using OpenFOAM Reader
reader = OpenFOAMReader("TES/TES.foam", cell_type=10)
velocity_func = reader.create_dolfinx_function(t=time, name="U")
foam_mesh = reader.dolfinx_mesh
#Get Tags from Mesh Regions
boundary = reader.reader.read()["boundary"]
inlet_facets, inlet_tags = tag_boundary_patch(foam_mesh, boundary["inlet"], 1)
outlet_facets, outlet_tags = tag_boundary_patch(foam_mesh, boundary["outlet"], 2)
vacuum_facets, vacuum_tags = tag_boundary_patch(foam_mesh, boundary["vacuum"], 4)
all_facets = np.concatenate([inlet_facets, outlet_facets, vacuum_facets])
all_tags = np.concatenate([inlet_tags, outlet_tags, vacuum_tags])
facet_tags = meshtags(foam_mesh, foam_mesh.topology.dim - 1, all_facets, all_tags)
# Create FESTIM mesh and assign region tags
mesh = F.Mesh(foam_mesh)
mesh._facet_tags = facet_tags
# Define Material (PbLi)
D_LiPb = htm.diffusivities.filter(material=htm.LIPB).mean()
material = F.Material(
D_0=D_LiPb.pre_exp.magnitude,
E_D=D_LiPb.act_energy.magnitude,
)
species = F.Species(material)
#Define Subdomains
fluid = F.VolumeSubdomain(id = 0, material = material)
inlet = F.SurfaceSubdomain(id = 1)
outlet = F.SurfaceSubdomain(id = 2)
vacuum = F.SurfaceSubdomain(id = 4)
volume_tags = mesh.define_meshtags(surface_subdomains = [], volume_subdomains=[fluid])
subdomains = [fluid, inlet, outlet, vacuum]
#Define Temperature - Function input
temperature = T
#Define Boundary Conditions
inlet_BC = F.FixedConcentrationBC(
subdomain=inlet,
value=1e15,
species=species
)
vacuum_BC = F.FixedConcentrationBC(
subdomain=vacuum,
value=0,
species=species
)
boundary_conditions = [inlet_BC, vacuum_BC]
#Define Simulation Settings
settings = F.Settings(
atol=1e-10,
rtol=1e-10,
transient=False,
)
#Define Advection Term
advection_term = [F.AdvectionTerm(velocity_func, fluid, species)]
#Define Desired Exports
c_field = F.XDMFExport(filename="H_concentration.xdmf", field=[species])
c_in = F.AverageSurface(species, inlet)
c_out = F.AverageSurface(species, outlet)
exports = [c_field, c_in, c_out]
my_sim = F.HydrogenTransportProblem(
mesh = mesh,
subdomains = subdomains,
species = [species],
temperature = temperature,
boundary_conditions = boundary_conditions,
settings = settings,
exports = exports,
advection_terms = advection_term
)
my_sim.initialise()
my_sim.ds = ufl.Measure("ds", domain=mesh.mesh, subdomain_data=facet_tags)
my_sim.run()
solver = my_sim.solver
return my_sim, c_in, c_out
my_sim, c_in, c_out = FESTIM(473.15, 1523)
c_out_value = c_out.value
c_in_value = c_in.value
print('c_in = ', c_in_value)
print('c_out = ', c_out_value)
#eta = 1 - (c_out_value/c_in_value)
#print('Extraction Efficiency = ', eta)