Diffusion of a Gaussian function - source term on for fraction of simulation time

dolfinx
1

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.

2

The issue is that:

  1. The magnitude of your source is very small, so it is not going to fundamentally change the solution.
  2. You are only applying the source for the first 0.1 seconds (5 time steps, so it doesn’t have much time to influence the solution).

Either:

  1. Increase the magnitude of the source
  2. Increase the time-span of the influence of the source.

Please note that you could do this with a single form, i.e.:

# From https://jsdokken.com/dolfinx-tutorial/chapter2/diffusion_code.html

import matplotlib as mpl
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.5

# 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 = uh.copy()
f.name = "uh"
f.x.array[:] *= 5.
L = (u_n + dt * f) * v * ufl.dx
bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = create_vector(linear_form)

solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)


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:
        f.x.array[:] = 0.
    assemble_vector(b, linear_form)

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

Awesome–that solves it. Thank you so much, and I appreciate the tip about using a single form!

4

Hello Fenics users,

I would like to slightly modify the tutorial Diffusion of a Gaussian function as
Gnordin did. In my case f = 0, u_{D} = 1 at the right and top boundaries. The initial condition here is u_{0} = 1.

In the code below, the changes from the tutorial code are on lines 33-35 for the initial condition (u (0, x, y) = 1) and on lines 59-67 for the boundary conditions. However, I don’t see the application of the boundary conditions (right and top boundary). What am I doing wrong? Thank you very much!

5

Here is my MWE:

# From https://jsdokken.com/dolfinx-tutorial/chapter2/diffusion_code.html

import matplotlib as mpl
import pyvista
import ufl
import numpy as np

from petsc4py import PETSc
from mpi4py import MPI

from dolfinx import fem, mesh, io, plot
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological, set_bc)
from dolfinx.io import XDMFFile, VTXWriter

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

# Define temporal parameters
t = 0  # Start time
T = 1.0  # Final time
num_steps = 50
dt = T / num_steps  # time step size

# Define mesh
nx, ny = 8, 8
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
                               [nx, ny], mesh.CellType.triangle)
V = fem.functionspace(domain, ("Lagrange", 1))

# Create initial condition u(0, x, y) = 1
def initial_condition(x, a=0):
    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)

#We create two python functions for determining the facets to apply boundary conditions to



def right(x):
    return np.isclose(x[0], 1)  #right (x = 1)


def top(x):
    return np.isclose(x[1], 1)  #top (y = 1)

# 
###################################################
######## Dirichlet boundary conditions #############
###################################################
# 

boundary_facets_right = locate_entities_boundary(domain, domain.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_right)
bc_m_right = dirichletbc(PETSc.ScalarType(0.), boundary_dofs_right, V)

boundary_facets_top = locate_entities_boundary(domain, domain.topology.dim - 1, top)
boundary_dofs_top = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_top)
bc_m_top = dirichletbc(PETSc.ScalarType(0.), boundary_dofs_top, V)

bcs = [bc_m_right, bc_m_top]  #dirichlet boundary conditions


#Time-dependent output
#xdmf = io.XDMFFile(domain.comm, "diffusion.xdmf", "w")
#xdmf.write_mesh(domain)

file = XDMFFile(domain.comm, "2DdiffusionEquation_june29th_2024.xdmf", "w")
file.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)
file.write_function(uh, t)


#Variational problem and solver
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))
a = u * v * ufl.dx + dt * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = (u_n + dt * f) * v * ufl.dx


#Preparing linear algebra structures for time dependent problems
bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = assemble_matrix(bilinear_form, bcs=bcs)
A.assemble()
b = create_vector(linear_form)


#Using petsc4py to create a linear solver
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)


#Visualization of time dependent problem using pyvista
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)])


#Updating the solution and right hand side per time step
for i in range(num_steps):
    t += dt

    # Update the right hand side reusing the initial vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    assemble_vector(b, linear_form)

    # Apply Dirichlet boundary condition to the vector
    apply_lifting(b, [bilinear_form], [bcs])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b, bcs)

    # 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
    file.write_function(uh, t)
    # Update plot
    new_warped = grid.warp_by_scalar("uh", factor=1)
    warped.points[:, :] = new_warped.points
    warped.point_data["uh"][:] = uh.x.array
    plotter.write_frame()
plotter.close()
file.close()
6

YOur bcs are not reflected in the code:

as both bcs are set to 0.

7

Thanks @dokken for your time.

Now let us modify the problem: f = 0, \ u_{D} = 1 at the right and top boundary. The initial condition is u(0, x, y) = 0.

Please tell me what this means?
Do I need to change the 0 value?

Thank you for your time and patience.

8

Here is my MWE :

# From https://jsdokken.com/dolfinx-tutorial/chapter2/diffusion_code.html

import matplotlib as mpl
import pyvista
import ufl
import numpy as np

from petsc4py import PETSc
from mpi4py import MPI

from dolfinx import fem, mesh, io, plot
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological, set_bc)
from dolfinx.io import XDMFFile, VTXWriter

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

# Define temporal parameters
t = 0  # Start time
T = 1.0  # Final time
num_steps = 100
dt = T / num_steps  # time step size

# Define mesh
nx, ny = 50, 50
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
                               [nx, ny], mesh.CellType.triangle)
V = fem.functionspace(domain, ("Lagrange", 1))

# Create initial condition u(0, x, y) = 1
def initial_condition(x, a=0):
    return np.exp(-a * (x[0]**2 + x[1]**2)) - 1


u_n = fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(initial_condition)

#We create two python functions for determining the facets to apply boundary conditions to



def right(x):
    return np.isclose(x[0], 1)  #right (x = 1)


def top(x):
    return np.isclose(x[1], 1)  #top (y = 1)

# 
###################################################
######## Dirichlet boundary conditions #############
###################################################
# 

boundary_facets_right = locate_entities_boundary(domain, domain.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_right)
bc_m_right = dirichletbc(PETSc.ScalarType(1.), boundary_dofs_right, V)

boundary_facets_top = locate_entities_boundary(domain, domain.topology.dim - 1, top)
boundary_dofs_top = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_top)
bc_m_top = dirichletbc(PETSc.ScalarType(1.), boundary_dofs_top, V)

bcs = [bc_m_right, bc_m_top]  #dirichlet boundary conditions


#Time-dependent output
#xdmf = io.XDMFFile(domain.comm, "diffusion.xdmf", "w")
#xdmf.write_mesh(domain)

file = XDMFFile(domain.comm, "2DdiffusionEquation_june30th_2024.xdmf", "w")
file.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)
file.write_function(uh, t)


#Variational problem and solver
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))
a = u * v * ufl.dx + dt * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = (u_n + dt * f) * v * ufl.dx


#Preparing linear algebra structures for time dependent problems
bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = assemble_matrix(bilinear_form, bcs=bcs)
A.assemble()
b = create_vector(linear_form)


#Using petsc4py to create a linear solver
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)


#Visualization of time dependent problem using pyvista
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)])


#Updating the solution and right hand side per time step
for i in range(num_steps):
    t += dt

    # Update the right hand side reusing the initial vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    assemble_vector(b, linear_form)

    # Apply Dirichlet boundary condition to the vector
    apply_lifting(b, [bilinear_form], [bcs])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b, bcs)

    # 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
    file.write_function(uh, t)
    # Update plot
    new_warped = grid.warp_by_scalar("uh", factor=1)
    warped.points[:, :] = new_warped.points
    warped.point_data["uh"][:] = uh.x.array
    plotter.write_frame()
plotter.close()
file.close()

but the results are not correct.