Error running steady-state advection only sim

User Support
1

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)
2

Welcome @charlieyates :tada:

   mesh._facet_tags = facet_tags

You should not interfere with “private” variables/attributes/methods (ie. starting with a _)

I think there are several problems here:

  • you don’t set the tags properly, the meshtags should be given to the problem class and not to the mesh
  • you shouldn’t have to use mesh.define_mesthags i believe

That would explain why your BC is not applied correctly.

You could try and export your meshtags to XDMF or VTX and then visualise them in Paraview, I believe you will see that your inlet surface isn’t tagged properly. I’m suspecting that the ordering of the mesh cells and facets in openFOAM is not the same as the one used in FESTIM, so you are tagging some cells but random ones in your geometry. Does that make sense?

tagging @jhdark

3

I’ve exported the dolfinx.mesh.Mesh with the meshtags to and XDMF and have viewed in ParaView where the tagging seems to be correct. Where should I give the meshtags in the problem if not to the mesh?

Below should be a screenshot of the tagged mesh in ParaView. Inlet should have id=1

tagged mesh
tagged mesh1264×728 104 KB

1 Like
4

Ok, so good! That means the tagging works.

I’ve been able to run your case (without advection for now) with:


import festim as F
import h_transport_materials as htm

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)

    # 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)

    # define meshtags but ignore surface subdomains
    _, 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_field =  F.VTXSpeciesExport(filename="H_concentration.bp", 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.facet_meshtags = facet_tags
    my_sim.volume_meshtags = volume_tags

    my_sim.initialise()
    my_sim.run()
    solver = my_sim.solver

    
    return my_sim, c_in, c_out

image
image2638×1251 94.9 KB

I’ll try to add the advection and see if that works

5

I’ve just tried with advection myself. No longer getting zero however c_out > c_in thus the calculated efficiency is negative.

6

Now I think the problem is a mesh issue.

Here is the h concentration field (with a capped colourbar for better visualisation) at 400 degC

image
image1920×1080 73.9 KB

You see that you have pretty high numerical inaccuracies next to the inlet.

That has to do with the local mesh Peclet number Pe∆ = U h/D

where U is the velocity, h is the cell size, and D is the diffusivity.

If you artificially increase the diffusivity by a factor 100, everything is much smoother:

image
image1920×1107 71.7 KB

And the efficiency is: 10%

7

Brilliant. Thank you for your help. I’m now getting 6.4%. How are you viewing the H_conc in ParaView? I’m exporting as XDMF but when I try open ParaView crashes.

8

I’m not exporting the XDMF but the VTX file

    c_field =  F.VTXSpeciesExport(filename="H_concentration.bp", field=[species])
9

@charlieyates just circling back on this, is it turbulent flow?