# Symmetry boundary condition **URL:** \ **Category:** User Support\ **Created:** [October 21, 2025, 9:45am UTC](https://festim.discourse.group/t/symmetry-boundary-condition/147 "2025-10-21T09:45:01Z")\ **Posts on this page:** 3\ **Page:** 1
**Author:** ![Andrea](https://avatars.discourse-cdn.com/v4/letter/a/c37758/32.png) [@Andrea](https://festim.discourse.group/u/Andrea)\ **Post date:** [October 21, 2025, 9:45am UTC](https://festim.discourse.group/t/symmetry-boundary-condition/147/1 "2025-10-21T09:45:01Z")
Hi all, I am working with a model that includes three different materials: tungsten, Eurofer, and FLiBe. The boundary conditions are shown in the picture below: ![image](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/45c4de6bb8b1e7a09c0fc1b19f0ea74890a8aa9b.jpeg) I am trying to simulate an advection–diffusion problem in which the FLiBe moves with a parabolic velocity profile. However, when I plot the mobile concentration, I notice that, instead of exiting through the top of the FLiBe region as expected, the concentration appears to leave through the Eurofer–FLiBe interface. ![image](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/29cd014cad37d3bf4fb1351afa72ab34858592ad.png) I suspect that this issue arises because no boundary condition has been imposed at the top of the FLiBe region. By default, the software applies a Neumann boundary condition with zero flux, which likely prevents the concentration from exiting there. To address this, I tried imposing a Neumann boundary condition using the known concentration in the FLiBe (previous plot 2.832e+19 [H/m^3]) and the average velocity (1.33 [m/s]). However, the resulting solution is affected by numerical errors: > Neum\_con=F.FluxBC(surfaces=12, value=2.832e+19\*1.33, field=0) #[H/m^2/s] ![image](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/88839ba359a4bc323caa8424b5f18716c298147b.jpeg) I am not sure how to correctly impose a symmetry boundary condition so that the concentration can leave from the top boundary. Could you please help me with this? I’m also attaching the code in case someone would like to reproduce the issue. Thank you in advance for your help. Best regards. [mesh\_domains.h5](https://festim.discourse.group/uploads/short-url/boyh5plMLgXUzpvB723yhEt2d4A.h5) (326.7 KB) [mesh\_domains.xdmf](https://festim.discourse.group/uploads/short-url/mMiu5uUu7x4ue4j9xFLEJ84bfA9.xdmf) (585 Bytes) [mesh\_boundaries.xdmf](https://festim.discourse.group/uploads/short-url/7cevLCEUZRXKwFe2fsiuEFNDCjC.xdmf) (588 Bytes) [mesh\_boundaries.h5](https://festim.discourse.group/uploads/short-url/mYV9sVDpIDzWmAvO0Xgas1vTuKi.h5) (195.1 KB) > import numpy as np > import meshio > import matplotlib.pyplot as plt > import festim as F > import os > import copy > from fenics import plot > > my\_model = F.Simulation() > > my\_model.mesh = F.MeshFromXDMF(volume\_file=“…/mesh\_domains.xdmf”, boundary\_file=“…/mesh\_boundaries.xdmf”) > > tungsten = F.Material( > id=8, > 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=9, > D\_0=4.57e-7/(np.sqrt(3)), > E\_D=0.23, > S\_0=1.35e22, > E\_S=0.16, > ) > > FLiBe = F.Material( > id=10, > 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] > from festim import x > import sympy as sp > > m1 = (1308.1 - 1337.0) / 0.01 # tra 0 e 0.01 > m2 = (1273.0 - 1308.1) / 0.03 # tra 0.01 e 0.04 > > T\_expr = sp.Piecewise( > (1337.0 + m1 \* x, x \<= 0.01), > (1308.1 + m2 \* (x - 0.01), (x \> 0.01) & (x \<= 0.04)), > (1273.0, True) > ) > > my\_model.T = F.Temperature(T\_expr) > > implatantion\_flux=F.ImplantationDirichlet(surfaces=6,phi=7.15e12, R\_p=9.52e-10, D\_0=tungsten.D\_0, E\_D=tungsten.E\_D) > > Neum\_con=F.FluxBC(surfaces=12, value=2.832e+19\*1.33, field=0) > > Diriclet\_con=F.DirichletBC(surfaces=7, value=4.6e-5\*(6.022e23), field=0) > > my\_model.boundary\_conditions=[implatantion\_flux, Diriclet\_con, Neum\_con] > > my\_model.sources=[F.Source(volume=10, field=0, value=1.1e19)] > > my\_model.settings=F.Settings( > absolute\_tolerance=1e-6, > relative\_tolerance=1e-6, > transient=False, > #final\_time=1000 > ) > “”" > my\_model.dt=F.Stepsize( > initial\_value=1e-4, > stepsize\_change\_ratio=1.2, > max\_stepsize=1, > dt\_min=1e-9, > milestones=[10, 100, 500, 1000] > ) > “”" > 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() > from fenics import interpolate, UserExpression, VectorFunctionSpace > > functionspace = VectorFunctionSpace(my\_model.mesh.mesh, “CG”, 1) > > class ParabolicVelocity(UserExpression): > def eval\_cell(self, value, x, ufc\_cell): > # x[0] = coordinata x, x[1] = y > if x[0] \< 0.04 or x[0] \> 0.06: > value[0] = 0.0 # componente x della velocità > value[1] = 0.0 # componente y (verticale) > else: > vmax = 2.0 > xc = 0.05 > r = 0.01 # raggio del profilo parabolico > vy = vmax \* (1 - ((x[0] - xc) / r) \*\* 2) > value[0] = 0.0 # flusso verticale, nessuna velocità orizzontale > value[1] = vy # componente verticale > def value\_shape(self): > return (2,) # perché è un vettore (vx, vy) > > velocity = ParabolicVelocity(degree=2) > > velocity = interpolate(velocity, functionspace) > CS = plot(velocity) > plt.colorbar(CS) > plt.show() > from dolfin import XDMFFile > > folder = “…” > > filepath = os.path.join(folder, “velocity\_field.xdmf”) > > xdmf = XDMFFile(filepath) > xdmf.write(velocity) > xdmf.close() > > from fenics import inner, dot, grad > my\_derived\_quantities = F.DerivedQuantities( > [ F.TotalVolume(field=“solute”, volume=10), > F.AverageVolume(field=“solute”, volume=10), > F.PointValue(field=“solute”, x=[0.05, 0.1]),]) > > my\_model.exports = [my\_derived\_quantities] > > my\_model.initialise() # reinitialisation is needed > > hydrogen\_concentration = my\_model.h\_transport\_problem.mobile.solution > test\_function\_solute = my\_model.h\_transport\_problem.mobile.test\_function > > advection\_term = inner(dot(grad(hydrogen\_concentration), velocity), test\_function\_solute) \* my\_model.mesh.dx(10) > > my\_model.h\_transport\_problem.F += advection\_term > print(my\_model.h\_transport\_problem.F) > > my\_model.run() > cp=plot(hydrogen\_concentration,title=" H transport mobile advection diffusion") > plt.colorbar(cp) > #plot(velocity) > plt.show() ---
**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:** [October 21, 2025, 2:38pm UTC](https://festim.discourse.group/t/symmetry-boundary-condition/147/2 "2025-10-21T14:38:33Z")
Cool model! In this problem, the right BC at the top/outlet of the FLiBe volume is a “do-nothing” BC so putting the default one is correct. An homogeneous Neumann BC enforces the gradient of concentration to be zero, that doesn’t mean that H isotopes cannot leave, the total flux is the diffusion flux + advective flux. The first screenshot you show makes sense to me. To convince yourself, you can try and calculate the inlet flux vs outlet advective fluxes with: J = \int\_{\Gamma} c \ v \cdot \mathbf{n} \ ds You can try and calculate with a back of the envelope calculation or by writing your own custom `DerivedQuantity` class similar to what is done in [this tutorial](https://festim-workshop.readthedocs.io/en/latest/content/task16.html#compute-the-advection-flux) ---
**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:** [October 21, 2025, 2:39pm UTC](https://festim.discourse.group/t/symmetry-boundary-condition/147/3 "2025-10-21T14:39:30Z")
Also, for this sort of application, I strongly recommend using FESTIM2 which comes with nice features for multi-domain and multi physics problems.