Time-stepping of nonlinear time dependent problem with time-dependent BC

Skonki · September 27, 2024, 2:00pm

Hello

I am trying to solve a heat-diffusion problem with a Robin boundary condition. For an evaluation at a given time I manage to solve the problem without any issue. However, as soon as I try to perform a time-stepping algorithm I get an infinite loop where the topolgy computation tries to find the connectivity:

2024-09-27 15:45:59.163 (1787.682s) [main ]topologycomputation.cpp:799 INFO| Requesting connectivity (2,0) - (1,0)
2024-09-27 15:45:59.163 (1787.683s) [main ]topologycomputation.cpp:799 INFO| Requesting connectivity (1,0) - (2,0)

The lines above repeat in an infinite loop.

I am using the latest release of dolfinx-fenics (https://fenicsproject.org) on a Lenovo Think Pad with an Intel(R) Core™ i7-8550U CPU @1.80GHZ. It’s a 64-Bit-System with Windows 11 Pro version 23H2. I have an Intel(R) UHD Graphics 620 graphic card.

And the code I am running is as follows and the error occures in the time-loop when I instatiate problem = NonlinearProblem(L, Tnew, bcs=bcs)

# Importing the functino space and numpy
import ufl
import numpy as np
import matplotlib.pyplot as plt
import pyvista
pyvista.set_jupyter_backend('html')

# Importing the solver bases
from mpi4py import MPI
from petsc4py import PETSc

# Importing the packages of fenicsx
from dolfinx import mesh, fem, io, log, default_scalar_type, plot
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

# Creating the mesh and the function space
domain = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100)
x = ufl.SpatialCoordinate(domain)                           # x is the 2D coordiante vector
V = fem.functionspace(domain, ("Lagrange", 1))

# Creating the trial function and the test function
T = fem.Function(V)                     # Function for Temperature (necessary as the problem is non-linear)
v = ufl.TestFunction(V)                 # Test function

# Defining numpy functions to find the boundaries
boundaries = [(1, lambda x: np.logical_or(np.isclose(x[0], 0),np.isclose(x[0], 1), np.isclose(x[1], 1))),
              (2, lambda x: np.isclose(x[1], 0))]

# Seperate as a function of the type
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = mesh.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 = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with io.XDMFFile(domain.comm, "facet_tags.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facet_tag, domain.geometry)

# Define the surface increment
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

# Defining the outward normal of the boundaries
n = ufl.FacetNormal(domain)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = fem.Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = fem.locate_dofs_topological(V, fdim, facets)
            self._bc = fem.dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = ufl.inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * ufl.inner(T ** 4 - values[1] ** 4, v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

# we define the Infinity distance air temperature function
T_inf_ex = lambda t: np.cos(2 * np.pi * (t - 17) / 24) * 10 + 7.5     # Daily variation [C°]

# define the boundary temperature
T_fix = lambda x: T_inf_ex(5) * x[0] ** 0 * x[1] ** 0                  # Boundary temp fix [C°]

# Defining the coefficients
epsilon = fem.Constant(domain, default_scalar_type(1))         # Thermal emmisivity      
sigma = fem.Constant(domain, default_scalar_type(100))         # Stefan-Boltzmann constant
k = fem.Constant(domain, default_scalar_type(1))               # Thermal conductivity
rho = fem.Constant(domain, default_scalar_type(1))             # Material density  
c = fem.Constant(domain, default_scalar_type(1))               # Specific heat capacity
Q = fem.Constant(domain, default_scalar_type(1E-5))            # Internal Heat generation

# Define the Dirichlet condition
boundary_conditions = [BoundaryCondition("Dirichlet", 1, T_fix),
                       BoundaryCondition("Robin", 2, (sigma * epsilon, T_inf_ex(0)))]

# Combining the system
L = k * ufl.inner(ufl.grad(T), ufl.grad(v)) * ufl.dx - ufl.inner(Q, v) * ufl.dx

# Adding Dirichlet bc
bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        L += condition.bc

# Solving the non-linearvariational problem
problem = NonlinearProblem(L, T, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
opts[f"{option_prefix}ksp_rtol"] = "1.e-8"
opts[f"{option_prefix}pc_type"] = "hypre"
opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
opts[f"{option_prefix}pc_hypre_boomeramg_max_iter"] = 1
opts[f"{option_prefix}pc_hypre_boomeramg_cycle_type"] = "v"
ksp.setFromOptions()

log.set_log_level(log.LogLevel.INFO)
niter, converged = solver.solve(T)
assert (converged)

# Define temporal parameters
t = to = 0  # Start time
tF = 96 # Final time
dt = 0.5
num_steps = int((tF - to) / 96)

xdmf = io.XDMFFile(domain.comm, "diffusion.xdmf", "w")
xdmf.write_mesh(domain)

# Define solution variable, and interpolate initial solution for visualization in Paraview
xdmf.write_function(T, to)

Tnew = T

for i in range(num_steps):
    t += dt
    Told = Tnew
    V = fem.functionspace(domain, ("Lagrange", 1))
    Tnew = fem.Function(V)
    v = ufl.TrialFunction(V)

    # Define the Dirichlet condition
    boundary_conditions = [BoundaryCondition("Dirichlet", 1, T_fix),
                           BoundaryCondition("Robin", 2, (dt * sigma * epsilon, T_inf_ex(t)))]

    # Combining the system
    L = dt * k * ufl.inner(ufl.grad(Tnew), ufl.grad(v)) * ufl.dx - dt * ufl.inner(Q, v) * ufl.dx  + rho * c * ufl.inner(Tnew, v) * ufl.dx - rho * c * np.dot(Told.x.array, v) * ufl.dx

    # Adding Dirichlet bc
    bcs = []
    for condition in boundary_conditions:
        if condition.type == "Dirichlet":
            bcs.append(condition.bc)
        else:
            L += condition.bc
    
    problem = NonlinearProblem(L, Tnew, bcs=bcs)
    solver = NewtonSolver(MPI.COMM_WORLD, problem)
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-6
    solver.report = True

    ksp = solver.krylov_solver
    opts = PETSc.Options()
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}ksp_type"] = "gmres"
    opts[f"{option_prefix}ksp_rtol"] = "1.e-8"
    opts[f"{option_prefix}pc_type"] = "hypre"
    opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
    opts[f"{option_prefix}pc_hypre_boomeramg_max_iter"] = 1
    opts[f"{option_prefix}pc_hypre_boomeramg_cycle_type"] = "v"
    ksp.setFromOptions()
    log.set_log_level(log.LogLevel.INFO)
    niter, converged = solver.solve(Tnew)
    assert (converged)
    xdmf.write_function(Tnew.x.array, to + (i+1) * dt)

xdmf.close()

Can anyone see where the problem lies?

francesco-ballarin · September 27, 2024, 2:56pm

One possible mistake is that in your BoundaryCondition class, type == "Robin" case, you seem to be using T, which is the unknown for the steady state solution. The unknown for the unsteady solution is called Tnew.

You could keep using T also for the unsteady case. Just define a Told = fem.Function(V) outside of the for i in range(num_steps) loop, and at every iteration copy whatever T or Tnew currently has into Told.

Skonki · September 28, 2024, 6:36am

Hello Francesco

Thank you very much for your suggestion.

I had already a similar setup previously with no luck. I have adapted the code as follows, the output remains the same.


Told = fem.Function(V)

for i in range(num_steps):
    t += dt
    Told = T
    
    T = fem.Function(V)
    v = ufl.TrialFunction(V)

    # Define the Dirichlet condition
    boundary_conditions = [BoundaryCondition("Dirichlet", 1, T_fix),
                           BoundaryCondition("Robin", 2, (dt * sigma * epsilon, T_inf_ex(t)))]

    # Combining the system
    L = dt * k * ufl.inner(ufl.grad(T), ufl.grad(v)) * ufl.dx - dt * ufl.inner(Q, v) * ufl.dx  + rho * c * ufl.inner(T, v) * ufl.dx - rho * c * np.dot(Told.x.array, v) * ufl.dx

    # Adding Dirichlet bc
    bcs = []
    for condition in boundary_conditions:
        if condition.type == "Dirichlet":
            bcs.append(condition.bc)
        else:
            L += condition.bc
    
    problem = NonlinearProblem(L, T, bcs=bcs)
    solver = NewtonSolver(MPI.COMM_WORLD, problem)
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-6
    solver.report = True

    ksp = solver.krylov_solver
    opts = PETSc.Options()
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}ksp_type"] = "gmres"
    opts[f"{option_prefix}ksp_rtol"] = "1.e-8"
    opts[f"{option_prefix}pc_type"] = "hypre"
    opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
    opts[f"{option_prefix}pc_hypre_boomeramg_max_iter"] = 1
    opts[f"{option_prefix}pc_hypre_boomeramg_cycle_type"] = "v"
    ksp.setFromOptions()
    log.set_log_level(log.LogLevel.INFO)
    niter, converged = solver.solve(Tnew)
    assert (converged)
    xdmf.write_function(T.x.array, to + (i+1) * dt)

xdmf.close()

francesco-ballarin · September 28, 2024, 6:06pm

Please avoid redefining your fem.Functions. Use one T and one Told, defined outside the loop, and update the content of the T.x.array.

Skonki · October 7, 2024, 12:43pm

Dear Francesco

this was what I had before. I have tried it again now and the switch does not change anything to the outcome. The problem remains.

francesco-ballarin · October 8, 2024, 5:16am

Please post the updated code with the above changes. The more you are able to simplify it, the higher the probabiltiy someone will be able to have a look in detail.

Skonki · October 8, 2024, 9:36am

Here the updated code as currently in use:

# Importing the functino space and numpy
import ufl
import numpy as np
import matplotlib.pyplot as plt
import pyvista
pyvista.set_jupyter_backend('html')

# Importing the solver bases
from mpi4py import MPI
from petsc4py import PETSc

# Importing the packages of fenicsx
from dolfinx import mesh, fem, io, log, default_scalar_type, plot
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

# Creating the mesh and the function space
domain = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100)
x = ufl.SpatialCoordinate(domain)                           # x is the 2D coordiante vector
V = fem.functionspace(domain, ("Lagrange", 1))

# Creating the trial function and the test function
T = fem.Function(V)                     # Function for Temperature (necessary as the problem is non-linear)
v = ufl.TestFunction(V)                 # Test function

# Defining numpy functions to find the boundaries
boundaries = [(1, lambda x: np.logical_or(np.isclose(x[0], 0),np.isclose(x[0], 1), np.isclose(x[1], 1))),
              (2, lambda x: np.isclose(x[1], 0))]

# Seperate as a function of the type
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = mesh.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 = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with io.XDMFFile(domain.comm, "facet_tags.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facet_tag, domain.geometry)

# Define the surface increment
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

# Defining the outward normal of the boundaries
n = ufl.FacetNormal(domain)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = fem.Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = fem.locate_dofs_topological(V, fdim, facets)
            self._bc = fem.dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = ufl.inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * ufl.inner(T ** 4 - values[1] ** 4, v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

# we define the Infinity distance air temperature function
T_inf_ex = lambda t: np.cos(2 * np.pi * (t - 17) / 24) * 10 + 7.5     # Daily variation [C°]

# define the boundary temperature
T_fix = lambda x: T_inf_ex(5) * x[0] ** 0 * x[1] ** 0                  # Boundary temp fix [C°]

# Defining the coefficients
epsilon = fem.Constant(domain, default_scalar_type(1))         # Thermal emmisivity      
sigma = fem.Constant(domain, default_scalar_type(100))         # Stefan-Boltzmann constant
k = fem.Constant(domain, default_scalar_type(1))               # Thermal conductivity
rho = fem.Constant(domain, default_scalar_type(1))             # Material density  
c = fem.Constant(domain, default_scalar_type(1))               # Specific heat capacity
Q = fem.Constant(domain, default_scalar_type(1E-5))            # Internal Heat generation

# Define the Dirichlet condition
boundary_conditions = [BoundaryCondition("Dirichlet", 1, T_fix),
                       BoundaryCondition("Robin", 2, (sigma * epsilon, T_inf_ex(0)))]

# Combining the system
L = k * ufl.inner(ufl.grad(T), ufl.grad(v)) * ufl.dx - ufl.inner(Q, v) * ufl.dx

# Adding Dirichlet bc
bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        L += condition.bc

# Solving the non-linearvariational problem
problem = NonlinearProblem(L, T, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
opts[f"{option_prefix}ksp_rtol"] = "1.e-8"
opts[f"{option_prefix}pc_type"] = "hypre"
opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
opts[f"{option_prefix}pc_hypre_boomeramg_max_iter"] = 1
opts[f"{option_prefix}pc_hypre_boomeramg_cycle_type"] = "v"
ksp.setFromOptions()

log.set_log_level(log.LogLevel.INFO)
niter, converged = solver.solve(T)
assert (converged)

# Define temporal parameters
t = to = 0  # Start time
tF = 96 # Final time
dt = 0.5
num_steps = int((tF - to) / 96)

xdmf = io.XDMFFile(domain.comm, "diffusion.xdmf", "w")
xdmf.write_mesh(domain)

# Define solution variable, and interpolate initial solution for visualization in Paraview
xdmf.write_function(T, to)

Told = fem.Function(V)

for i in range(num_steps):
    t += dt
    Told = T

    # Define the Dirichlet condition
    boundary_conditions = [BoundaryCondition("Dirichlet", 1, T_fix),
                           BoundaryCondition("Robin", 2, (dt * sigma * epsilon, T_inf_ex(t)))]

    # Combining the system
    L = dt * k * ufl.inner(ufl.grad(T), ufl.grad(v)) * ufl.dx - dt * ufl.inner(Q, v) * ufl.dx  + rho * c * ufl.inner(T, v) * ufl.dx - rho * c * np.dot(Told.x.array, v) * ufl.dx

    # Adding Dirichlet bc
    bcs = []
    for condition in boundary_conditions:
        if condition.type == "Dirichlet":
            bcs.append(condition.bc)
        else:
            L += condition.bc
    
    problem = NonlinearProblem(L, T, bcs=bcs)
    solver = NewtonSolver(MPI.COMM_WORLD, problem)
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-6
    solver.report = True

    ksp = solver.krylov_solver
    opts = PETSc.Options()
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}ksp_type"] = "gmres"
    opts[f"{option_prefix}ksp_rtol"] = "1.e-8"
    opts[f"{option_prefix}pc_type"] = "hypre"
    opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
    opts[f"{option_prefix}pc_hypre_boomeramg_max_iter"] = 1
    opts[f"{option_prefix}pc_hypre_boomeramg_cycle_type"] = "v"
    ksp.setFromOptions()
    log.set_log_level(log.LogLevel.INFO)
    niter, converged = solver.solve(Tnew)
    assert (converged)
    xdmf.write_function(T.x.array, to + (i+1) * dt)

xdmf.close()

francesco-ballarin · October 8, 2024, 11:01am

Fixed version that seems to work for me. Look for lines with # CHANGED comment.

# Importing the functino space and numpy
import ufl
import numpy as np

# Importing the solver bases
from mpi4py import MPI
from petsc4py import PETSc

# Importing the packages of fenicsx
from dolfinx import mesh, fem, io, log, default_scalar_type, plot
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

# Creating the mesh and the function space
domain = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100)
x = ufl.SpatialCoordinate(domain)                           # x is the 2D coordiante vector
V = fem.functionspace(domain, ("Lagrange", 1))

# Creating the trial function and the test function
T = fem.Function(V)                     # Function for Temperature (necessary as the problem is non-linear)
v = ufl.TestFunction(V)                 # Test function

# Defining numpy functions to find the boundaries
boundaries = [(1, lambda x: np.logical_or(np.isclose(x[0], 0),np.isclose(x[0], 1), np.isclose(x[1], 1))),
              (2, lambda x: np.isclose(x[1], 0))]

# Seperate as a function of the type
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = mesh.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 = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with io.XDMFFile(domain.comm, "facet_tags.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facet_tag, domain.geometry)

# Define the surface increment
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

# Defining the outward normal of the boundaries
n = ufl.FacetNormal(domain)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = fem.Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = fem.locate_dofs_topological(V, fdim, facets)
            self._bc = fem.dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = ufl.inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * ufl.inner(T ** 4 - values[1] ** 4, v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

# we define the Infinity distance air temperature function
T_inf_ex = lambda t: np.cos(2 * np.pi * (t - 17) / 24) * 10 + 7.5     # Daily variation [C°]

# define the boundary temperature
T_fix = lambda x: T_inf_ex(5) * x[0] ** 0 * x[1] ** 0                  # Boundary temp fix [C°]

# Defining the coefficients
epsilon = fem.Constant(domain, default_scalar_type(1))         # Thermal emmisivity
sigma = fem.Constant(domain, default_scalar_type(100))         # Stefan-Boltzmann constant
k = fem.Constant(domain, default_scalar_type(1))               # Thermal conductivity
rho = fem.Constant(domain, default_scalar_type(1))             # Material density
c = fem.Constant(domain, default_scalar_type(1))               # Specific heat capacity
Q = fem.Constant(domain, default_scalar_type(1E-5))            # Internal Heat generation

# Define the Dirichlet condition
boundary_conditions = [BoundaryCondition("Dirichlet", 1, T_fix),
                       BoundaryCondition("Robin", 2, (sigma * epsilon, T_inf_ex(0)))]

# Combining the system
L = k * ufl.inner(ufl.grad(T), ufl.grad(v)) * ufl.dx - ufl.inner(Q, v) * ufl.dx

# Adding Dirichlet bc
bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        L += condition.bc

# Solving the non-linearvariational problem
problem = NonlinearProblem(L, T, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
opts[f"{option_prefix}ksp_rtol"] = "1.e-8"
opts[f"{option_prefix}pc_type"] = "hypre"
opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
opts[f"{option_prefix}pc_hypre_boomeramg_max_iter"] = 1
opts[f"{option_prefix}pc_hypre_boomeramg_cycle_type"] = "v"
ksp.setFromOptions()

log.set_log_level(log.LogLevel.INFO)
niter, converged = solver.solve(T)
assert (converged)

# Define temporal parameters
t = to = 0  # Start time
tF = 96 # Final time
dt = 0.5
num_steps = int((tF - to) / 96)

xdmf = io.XDMFFile(domain.comm, "diffusion.xdmf", "w")
xdmf.write_mesh(domain)

# Define solution variable, and interpolate initial solution for visualization in Paraview
xdmf.write_function(T, to)

Told = fem.Function(V)

for i in range(num_steps):
    t += dt
    Told.x.array[:] = T.x.array  # CHANGED

    # Define the Dirichlet condition
    boundary_conditions = [BoundaryCondition("Dirichlet", 1, T_fix),
                           BoundaryCondition("Robin", 2, (dt * sigma * epsilon, T_inf_ex(t)))]

    # Combining the system
    L = dt * k * ufl.inner(ufl.grad(T), ufl.grad(v)) * ufl.dx - dt * ufl.inner(Q, v) * ufl.dx  + rho * c * ufl.inner(T, v) * ufl.dx - rho * c * ufl.inner(Told, v) * ufl.dx  # CHANGED

    # Adding Dirichlet bc
    bcs = []
    for condition in boundary_conditions:
        if condition.type == "Dirichlet":
            bcs.append(condition.bc)
        else:
            L += condition.bc

    problem = NonlinearProblem(L, T, bcs=bcs)
    solver = NewtonSolver(MPI.COMM_WORLD, problem)
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-6
    solver.report = True

    ksp = solver.krylov_solver
    opts = PETSc.Options()
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}ksp_type"] = "gmres"
    opts[f"{option_prefix}ksp_rtol"] = "1.e-8"
    opts[f"{option_prefix}pc_type"] = "hypre"
    opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
    opts[f"{option_prefix}pc_hypre_boomeramg_max_iter"] = 1
    opts[f"{option_prefix}pc_hypre_boomeramg_cycle_type"] = "v"
    ksp.setFromOptions()
    log.set_log_level(log.LogLevel.INFO)
    niter, converged = solver.solve(T)  # CHANGED
    assert (converged)
    xdmf.write_function(T, to + (i+1) * dt)  # CHANGED

xdmf.close()

Skonki · October 8, 2024, 12:31pm

Hello Francesco

this indeed solved the issue, now I have a convergence problem but that’s on my model setup side. Thank you very much.

Could you tell me why the update of the values of T (through T.x.array) instead of the entire function did the trick? It is a bit strange to me how this worked and I would like to learn from this mistake.

francesco-ballarin · October 8, 2024, 3:08pm

Could you tell me why the update of the values of T (through T.x.array) instead of the entire function did the trick

In ufl (and also in your BoundaryCondition class) you define forms using a dolfinx.fem.Function. Once T is in the ufl form (or in the BC), if you suddenly rename the python variable T to be someting else, the new variable will not replace the old one in forms/expressions where you had previously used the old variable.

SYSYSYSY · November 21, 2024, 7:37am

num_steps = int((tf - to )/dt)

Topic		Replies	Views
Time-dependent mesh and nonlinear problem (PETSC Error after migrating from FEniCSX 0.7.2 to 0.8.0) General	2	45	October 20, 2024
Alternative formulation for Time dependent Dirichlet BC dolfinx	16	358	January 2, 2024
Time-dependent nonlinear reaction-diffusion not proceeding in time loop	1	559	January 13, 2023
Error in implementing 2D time dependent heat equation dolfinx	4	172	March 6, 2024
Getting same results at timesteps 0 and 1 dolfinx	9	46	December 10, 2024

Time-stepping of nonlinear time dependent problem with time-dependent BC

Related topics