I want to make a conceptually simple modification of the tutorial Diffusion of a Gaussian function where the Gaussian initial condition is carried over as a source term for a fraction of the simulation time, say 0 \le t \le 0.1. In the code below, the changes from the tutorial code are on lines 64-69, 72-73, and 124-127. However, this code gives the same result as the unmodified tutorial code. What am I doing wrong? Thanks!
# From https://jsdokken.com/dolfinx-tutorial/chapter2/diffusion_code.html
import numpy as np
import pyvista
import ufl
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import (
apply_lifting,
assemble_matrix,
assemble_vector,
create_vector,
set_bc,
)
from mpi4py import MPI
from petsc4py import PETSc
# Define temporal parameters
t = 0 # Start time
T = 1.0 # Final time
num_steps = 50
dt = T / num_steps # time step size
t_on = 0.1
# Define mesh
nx, ny = 50, 50
domain = mesh.create_rectangle(
MPI.COMM_WORLD,
[np.array([-2, -2]), np.array([2, 2])],
[nx, ny],
mesh.CellType.triangle,
)
V = fem.FunctionSpace(domain, ("Lagrange", 1))
# Create initial condition
def initial_condition(x, a=5):
return np.exp(-a * (x[0] ** 2 + x[1] ** 2))
u_n = fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(initial_condition)
# Create boundary condition
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(
domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool)
)
bc = fem.dirichletbc(
PETSc.ScalarType(0), fem.locate_dofs_topological(V, fdim, boundary_facets), V
)
xdmf = io.XDMFFile(domain.comm, "diffusion.xdmf", "w")
xdmf.write_mesh(domain)
# Define solution variable, and interpolate initial solution for visualization in Paraview
uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)
xdmf.write_function(uh, t)
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
a = u * v * ufl.dx + dt * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
# Source off
f_source_off = fem.Constant(domain, PETSc.ScalarType(0))
L_source_off = (u_n + dt * f_source_off) * v * ufl.dx
# Source off
f_source_on = uh.copy()
L_source_on = (u_n + dt * f_source_on) * v * ufl.dx
bilinear_form = fem.form(a)
linear_form_source_off = fem.form(L_source_off)
linear_form_source_on = fem.form(L_source_on)
A = assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = create_vector(linear_form_source_off)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
import matplotlib as mpl
pyvista.start_xvfb()
grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(V))
plotter = pyvista.Plotter()
plotter.open_gif("u_time.gif", fps=10)
grid.point_data["uh"] = uh.x.array
warped = grid.warp_by_scalar("uh", factor=1)
viridis = mpl.colormaps.get_cmap("viridis").resampled(25)
sargs = dict(
title_font_size=25,
label_font_size=20,
fmt="%.2e",
color="black",
position_x=0.1,
position_y=0.8,
width=0.8,
height=0.1,
)
renderer = plotter.add_mesh(
warped,
show_edges=True,
lighting=False,
cmap=viridis,
scalar_bar_args=sargs,
clim=[0, max(uh.x.array)],
)
for i in range(num_steps):
t += dt
# Update the right hand side reusing the initial vector and appropriate linear
# form for source on or off
with b.localForm() as loc_b:
loc_b.set(0)
if t <= t_on:
assemble_vector(b, linear_form_source_on)
else:
assemble_vector(b, linear_form_source_off)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bc])
# Solve linear problem
solver.solve(b, uh.vector)
uh.x.scatter_forward()
# Update solution at previous time step (u_n)
u_n.x.array[:] = uh.x.array
# Write solution to file
xdmf.write_function(uh, t)
# Update plot
warped = grid.warp_by_scalar("uh", factor=1)
plotter.update_coordinates(warped.points.copy(), render=False)
plotter.update_scalars(uh.x.array, render=False)
plotter.write_frame()
plotter.close()
xdmf.close()
P.S. I have gone through related topics in the forum, but did not find information that helped me resolve the problem.