How to implement inhomogeneous Neumann boundary condition in FeniCSX?

dolfinx
1

Hi dear FeniCSX users,

I am using FeniCSX to solve the following 1D coupled system:

\partial_{\tilde{t}} \tilde{m}= \mathrm{Lu} \partial_{xx} \tilde{m} - \mathrm{Lu}\,\mathrm{Pn} \partial_{xx} \tilde{T}
\partial_{\tilde{t}} \tilde{T}= \partial_{xx} \tilde{T} - \varepsilon\mathrm{Ko} \partial_{\tilde{t}} \tilde{m}

with the boundary conditions

\partial_{\tilde{x}} \tilde{m} (0, \tilde{t}) - Pn \partial_{\tilde{x}} \tilde{T} (0, \tilde{t}) = 0
\partial_{\tilde{x}} \tilde{T} (0, \tilde{t}) = -Q

- \partial_{\tilde{x}} \tilde{m} ( 1, \tilde{t}) + Pn \partial_{\tilde{x}} \tilde{T} (1, \tilde{t}) + B_{im} (1 - \tilde{m}) = 0

\partial_{\tilde{x}} \tilde{T} ( 1, \tilde{t})= B_{it} \bigl(1 - \tilde{T}(1, \tilde{t}) \bigr) + (1 - \varepsilon) Ko Lu B_{im} [ \tilde{m} (1, \tilde{t}) - 1]

and the initial conditions

\tilde{m}(\tilde{x}, ,0)=0
\tilde{T}(\tilde{x} ,0)=0

The weak formulation of this problem is given by:

F0 = \int \partial_{\tilde{t}} \tilde{m} q + Lu \int \partial_{\tilde{x}} \tilde{m} \partial_{\tilde{x}} q - Lu Pn \int \partial_{\tilde{x}} \tilde{T} \partial_{\tilde{x}} q + Lu \int_{x=1} B_{im} [m-1] q

F1 = \int \partial_{\tilde{t}} \tilde{T} v + \int \partial_{\tilde{x}} \tilde{T} \partial_{\tilde{x}} v - \varepsilon Ko \int \partial_{\tilde{t}} \tilde{m} v + \int_{x=0} Q v + \int_{x=1} B_{it} \bigl( \tilde{T}(1, \tilde{t}) - 1 \bigr) - (1 - \varepsilon) Ko Lu B_{im} [ \tilde{m} (1, \tilde{t}) - 1] v

Below is my MWE with FeniCSX:
First attempt:


F0 =  inner(m, q) * dx - inner(m_n, q) * dx + dt * Lu_true * inner(grad(m), grad(q)) * dx - dt * ( Lu_true * Pn_true ) * inner(grad(T), grad(q)) * dx

Term_BCs_m = Lu_true * dt * inner(Bim_true*(m-1), q) * ds(2)    

F0_ =  F0 + Term_BCs_m

F1 = inner(T, v) * dx - inner(T_n, v) * dx + dt * inner(grad(T), grad(v)) * dx + (epsilon_true * Ko_true) * inner(m, v) * dx - (epsilon_true * Ko_true) * inner(m_n, v) * dx 

Term_BCs_T_c =  dt * inner(Q, v) * ds(1)  
Term_BCs_T_d =  dt * inner(Bit_true*(T-1) - (1-epsilon_true)*Ko_true*Lu_true*Bim_true(m-1), v) * ds(2)

F1_ = F1 + Term_BCs_T_c + Term_BCs_T_d

Second attempt:

n = FacetNormal(domain)
Q = -dot(n, grad(T))

F0 =  inner(m, q) * dx - inner(m_n, q) * dx + dt * Lu_true * inner(grad(m), grad(q)) * dx - dt * ( Lu_true * Pn_true ) * inner(grad(T), grad(q)) * dx

Term_BCs_m = Lu_true * dt * inner(Bim_true*(m-1), q) * ds(2)    

F0_ =  F0 + Term_BCs_m

F1 = inner(T, v) * dx - inner(T_n, v) * dx + dt * inner(grad(T), grad(v)) * dx + (epsilon_true * Ko_true) * inner(m, v) * dx - (epsilon_true * Ko_true) * inner(m_n, v) * dx 

Term_BCs_T_c =  dt * inner(-Q, v) * ds(1)  
Term_BCs_T_d =  dt * inner(Bit_true*(T-1) - (1-epsilon_true)*Ko_true*Lu_true*Bim_true(m-1), v) * ds(2)

F1_ = F1 + Term_BCs_T_c + Term_BCs_T_d

Both gave the same result. But the temperature should be between [-4, 2] but in my exceeds 2 and equal to 13.

I am struggling to understand how to define the inhomogeneous Neumann boundary condition in my case.

Any help would be appreciate.

Thank you.

2 

import numpy as np

import ufl
import basix

import dolfinx.fem.petsc


from dolfinx.fem import (dirichletbc, Function, FunctionSpace,
                          locate_dofs_topological)

from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities, meshtags, locate_entities_boundary
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, TrialFunction, derivative, Measure, SpatialCoordinate, FacetNormal, dot
from dolfinx import fem, io, mesh
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType

from dolfinx import log, default_scalar_type

from mpi4py import MPI

import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")



# parameters
xmin = 0  #the coordinate of left boundary
xmax = 1  #the coordinate of right boundary

# Define mesh            #nums =  [8, 16, 32,64,96, 128,255] number of cells
nx = 20 #100         # Modified on Tue 3 Feb 2026  
domain = mesh.create_interval(MPI.COMM_WORLD, nx, [xmin, xmax])      #Added on Mon 2 Feb 2026

x = SpatialCoordinate(domain)  


###################################################
#### Parameters : time and load and material ######
###################################################
# Case study : Ceramics


#Non-dimensional parameters

Ko_true = fem.Constant(domain, ScalarType(49.42))             
Pn_true = fem.Constant(domain, ScalarType(0.084))          
Lu_true = fem.Constant(domain, ScalarType(0.238))           
epsilon_true = fem.Constant(domain, ScalarType(0.2))        
Bim_true = fem.Constant(domain, ScalarType(3.33))              
Bit_true = fem.Constant(domain, ScalarType(2.5))             
#Qe = fem.Constant(domain, ScalarType(0.9))     
#Le_true = 1/Lu_true                


#Modified on Sat 22 Feb 2025

print('Lu_true =', Lu_true, 'Pn_true =', Pn_true, 'epsilontrue =', epsilon_true, 'Ko_true =', Ko_true, 'Bim_true=', Bim_true, 'Bit_true=', Bit_true )


#
###################################################
########## Mixed Space of resolution ##############
###################################################
# 

# Define ME function space
# Lagrange" is type of the finite element, and 1(piecewise linear elements) or 2 (piecewise quadratic elements) is the degree
m1 = basix.ufl.element("Lagrange", domain.basix_cell(), 1)             # For moisture distribution  # Modified on Wed 18 June 2025
T1 = basix.ufl.element("Lagrange", domain.basix_cell(), 1)             # For temperature evolution  # Modified on Wed 18 June 2025

mixed_el = basix.ufl.mixed_element([m1, T1])  # Mixed element

#Define the function space
ME = fem.functionspace(domain, mixed_el)     # Mixed Function space for all dofs


# Create the initial and solution functions of space
u = Function(ME)  # current solution  # total dofs in the problem/ func for all dofs
print("u dofs:", len(u.x.array))

m_sp = fem.functionspace(domain, m1) # moisture function space
T_sp = fem.functionspace(domain, T1) # temperature function space

Vm, uu_to_m = ME.sub(0).collapse() # dofs in u func subspace
m_n = fem.Function(Vm) # dofs in u func

VT, u_to_T = ME.sub(1).collapse() # dofs in u func subspace
T_n = fem.Function(VT) # dofs in u func

m, T = ufl.split(u) # (trail function) give two subfunctions for two elements/dofs from mixed space. Otherwise decoupled. for nonlinear don't use trial fn
(q, v) = ufl.TestFunctions(ME)


## # Define temporal parameters / Time parametrization
t         = 0.               # Start time
Tf        = 10.           # End time
num_steps = 1000    # Number of time steps 
dt = ( (Tf-t)/num_steps ) 



#----------------------Initial Conditions---------------------------#
# Initialize history values
u.x.array[:] = 0.0

m_n.x.array[:] = 0.0
T_n.x.array[:] = 0.0

#dimensionless initial moisture content
def m_init(x):
    values = np.zeros((1, x.shape[1]))
    values[0] = 0.0
    return values

#dimensionless initial temperature content
def T_init(x):
    values = np.zeros((1, x.shape[1]))
    values[0] = 0.0 
    return values

m_n.interpolate(m_init)
T_n.interpolate(T_init)

# checking if it is properly imlemented
m_n.name = "moisture"           #initial moisture
T_n.name = "temperature"        #initial temperature

 
#----------------------Boundary Conditions---------------------------#
# We start by identifying the facets contained in each boundary and create a
# custom integration measure ds
boundaries = [(1, lambda x: np.isclose(x[0], xmin)), # left (x=0)
              (2, lambda x: np.isclose(x[0], xmax))] # right (x=1)

# We now loop through all the boundary conditions and create MeshTags
# identifying the facets for each boundary condition.
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

# Debugging boundary condition
# To debug boundary conditions, the easiest thing to do is to visualize the boundary in Paraview
# by writing the MeshTags to file. We can then inspect individual boundaries using the Threshold-filter.
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with XDMFFile(domain.comm, "facet_tags_Thu9Apr2026_1D_Backward_Euler_Low_Order_Basis_Functions_ConvectiveBCs.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facet_tag, domain.geometry)



# Integration measures
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata) #line/area element
dx = ufl.Measure("dx", domain=domain, metadata=metadata)                           #volume element

 

Qe = fem.Constant(domain, ScalarType(0.9))       #prescribed heat flux
#n = FacetNormal(domain)
#Q = -dot(n, grad(T))
#print(Q)
F0 =  inner(m, q) * dx - inner(m_n, q) * dx + dt * Lu_true * inner(grad(m), grad(q)) * dx - dt * ( Lu_true * Pn_true ) * inner(grad(T), grad(q)) * dx

Term_BCs_m = Lu_true * dt * inner(Bim_true*(m-1), q) * ds(2)  


F0_ =  F0 + Term_BCs_m

F1 = inner(T, v) * dx - inner(T_n, v) * dx + dt * inner(grad(T), grad(v)) * dx + (epsilon_true * Ko_true) * inner(m, v) * dx - (epsilon_true * Ko_true) * inner(m_n, v) * dx

Term_BCs_T_c =  dt * inner(Qe, v) * ds(1)  
Term_BCs_T_d =  dt * inner(Bit_true*(T-1) - (1-epsilon_true)*Ko_true*Lu_true*Bim_true(m-1), v) * ds(2)

F1_ = F1 + Term_BCs_T_c + Term_BCs_T_d


F = F0_ + F1_


# Non linear problem definition
problem = NonlinearProblem(F, u)



#Solve variational problem
# Create nonlinear problem and Newton solver

  
solver = NewtonSolver(domain.comm, problem)                          #solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"      #residual    #incremental
# Absolute tolerance
solver.atol = 5e-10       # 1e-8
solver.rtol = 1e-12 #0.0 #1e-14 #0.0 #1e-12    # 1e-8
#solver.max_it = 70 
solver.report = True
#solver.max_it = 100        #To increase the number of iterations
solver.error_on_nonconvergence = False
solver.max_it = 3


#-----------------------------Solution Settings-----------------------------#
Vm_sol, u_to_m_sol = ME.sub(0).collapse() # dofs in u func subspace
m_sol = fem.Function(Vm_sol)

VT_sol, u_to_T_sol = ME.sub(1).collapse() # dofs in u func subspace
T_sol = fem.Function(VT_sol)

#m_sol.name = "moisture"
#T_sol.name = "temperature"

file_moist = dolfinx.io.VTXWriter(domain.comm, "1D_moisture_Backward_Euler_Low_Order_Basis_Functions_Thu9Apr2026_ConvectiveBCs_100mesh.bp", m_n, "bp4")
file_temp = dolfinx.io.VTXWriter(domain.comm, "1D_temperature_Backward_Euler_Low_Order_Basis_Functions_Thu9Apr2026_ConvectiveBCs.bp", T_n, "bp4")


import time

start_time = time.time()


#---------------------------Solve problem-------------------------------------#
duration_solve = 0.0
t = 0.0
file_moist.write(t)           # Create file for initial moisture
file_temp.write(t)            # Create file for initial temperature
# print 
num_steps = int(num_steps)
for n in range(num_steps):
    # Update time
    t += dt

    # Solve newton-steps
    log.set_log_level(log.LogLevel.INFO)
    duration_solve -= MPI.Wtime()
    niter, converged = solver.solve(u)
    duration_solve += MPI.Wtime()

    PETSc.Sys.Print(
        "Phys. Time {:.4f}, Calc. Time {:.4f}, Num. Iter. {}".format(
            t, duration_solve, niter
        )
    )

    m_n.x.array[:] = m_sol.x.array
    T_n.x.array[:] = T_sol.x.array

    m_sl = u.sub(0).collapse()
    T_sl = u.sub(1).collapse()

    m_sl.x.scatter_forward()
    T_sl.x.scatter_forward()

    m_sol.interpolate(m_sl)
    T_sol.interpolate(T_sl)

    m_n.x.array[:] = m_sol.x.array
    T_n.x.array[:] = T_sol.x.array

    file_moist.write(t)
    file_temp.write(t)

    

file_moist.close()
file_temp.close()

# outfile.close()


end_time = time.time()

print('final time is = ', end_time - start_time)

print(Q)

This is the code I wrote to solve the above system. However, I did not get the expected results. Any help would be much appreciated.