# Error running steady-state advection only sim

**URL:** <https://festim.discourse.group/t/error-running-steady-state-advection-only-sim/115>\
**Category:** User Support\
**Tags:** festim2\
**Created:** [April 24, 2025, 1:01pm UTC](https://festim.discourse.group/t/error-running-steady-state-advection-only-sim/115 "2025-04-24T13:01:19Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![charlieyates](https://avatars.discourse-cdn.com/v4/letter/c/71c47a/32.png) [@charlieyates](https://festim.discourse.group/u/charlieyates)\
**Post date:** [April 24, 2025, 1:01pm UTC](https://festim.discourse.group/t/error-running-steady-state-advection-only-sim/115/1 "2025-04-24T13:01:19Z")

</div>

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:

```auto
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

```

```auto

  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)

```

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [April 24, 2025, 2:51pm UTC](https://festim.discourse.group/t/error-running-steady-state-advection-only-sim/115/2 "2025-04-24T14:51:07Z")

</div>

Welcome @charlieyates 🎉

```python
   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

---

<div class="post-metadata">

**Author:** ![charlieyates](https://avatars.discourse-cdn.com/v4/letter/c/71c47a/32.png) [@charlieyates](https://festim.discourse.group/u/charlieyates)\
**Post date:** [April 24, 2025, 3:39pm UTC](https://festim.discourse.group/t/error-running-steady-state-advection-only-sim/115/3 "2025-04-24T15:39:56Z")

</div>

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](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/7f7554279ef7fcc11f19ca0e7511f01a23dff281.png)

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [April 24, 2025, 3:54pm UTC](https://festim.discourse.group/t/error-running-steady-state-advection-only-sim/115/4 "2025-04-24T15:54:50Z")

</div>

Ok, so good! That means the tagging works.

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

```python

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](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/06f609b4baa07cddcec26720d9b20658c796cd65.jpeg)

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

---

<div class="post-metadata">

**Author:** ![charlieyates](https://avatars.discourse-cdn.com/v4/letter/c/71c47a/32.png) [@charlieyates](https://festim.discourse.group/u/charlieyates)\
**Post date:** [April 24, 2025, 4:19pm UTC](https://festim.discourse.group/t/error-running-steady-state-advection-only-sim/115/5 "2025-04-24T16:19:51Z")

</div>

I’ve just tried with advection myself. No longer getting zero however c\_out \> c\_in thus the calculated efficiency is negative.

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [April 24, 2025, 4:30pm UTC](https://festim.discourse.group/t/error-running-steady-state-advection-only-sim/115/6 "2025-04-24T16:30:30Z")

</div>

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](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/eaf45e69471305785f69811143134f4f8aa94f90.jpeg)

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](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/d2928570727725a2be1de4333fbc234a5cb9386e.jpeg)

And the efficiency is: 10%

---

<div class="post-metadata">

**Author:** ![charlieyates](https://avatars.discourse-cdn.com/v4/letter/c/71c47a/32.png) [@charlieyates](https://festim.discourse.group/u/charlieyates)\
**Post date:** [April 24, 2025, 4:47pm UTC](https://festim.discourse.group/t/error-running-steady-state-advection-only-sim/115/7 "2025-04-24T16:47:49Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [April 24, 2025, 5:00pm UTC](https://festim.discourse.group/t/error-running-steady-state-advection-only-sim/115/8 "2025-04-24T17:00:03Z")

</div>

I’m not exporting the XDMF but the VTX file

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

```

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [June 13, 2025, 7:40pm UTC](https://festim.discourse.group/t/error-running-steady-state-advection-only-sim/115/9 "2025-06-13T19:40:29Z")

</div>

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