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