Time dependent diffusion constant in dolfinx

igarashi.takahiro · April 25, 2022, 4:38am

Hello everyone.

We would like to perform a diffusion calculation in which the diffusion coefficient changes with time.
We succeeded with the program shown below.
However, solver object is re-defined for each time step.
We think that is not efficient.
Can we conduct a calculation with varying diffusion coefficients without redefining the solver object?

# Module
from mpi4py import MPI
from dolfinx.mesh import create_rectangle
from dolfinx.fem import FunctionSpace, Function, Constant, locate_dofs_geometrical, dirichletbc, form, Expression
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import XDMFFile
from ufl import TrialFunction, TestFunction, dot, grad, dx, lhs, rhs
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import LinearProblem
import numpy as np

# Parameters
lx,ly,nx,ny = 10.0,5.0,10,5
vl,vr,vt,vb = 5.0,10.0,10.0,0.0

D = 1.0
f = 0.0
dt = 1.0
u0 = 0.0

# Mesh
domain = create_rectangle(MPI.COMM_WORLD, ([0, 0],[lx, ly]), [nx, ny])

# Function space
V = FunctionSpace(domain, ("CG", 1))

# Functions
u,v = TrialFunction(V),TestFunction(V)
u_n = Function(V)

# initial value
u_n.vector.array[:] = u0 * np.ones(u_n.vector.array.shape)

# Boundary condition
def left_boundary(x):
    return(np.isclose(x[0], 0.0))

def right_boundary(x):
    return(np.isclose(x[0], lx))

def top_boundary(x):
    return(np.isclose(x[1], ly))

def bottom_boundary(x):
    return(np.isclose(x[1], 0.0))

b_dofs_l = locate_dofs_geometrical(V, left_boundary)
b_dofs_r = locate_dofs_geometrical(V, right_boundary)
b_dofs_t = locate_dofs_geometrical(V, top_boundary)
b_dofs_b = locate_dofs_geometrical(V, bottom_boundary)

bc_l = dirichletbc(ScalarType(vl), b_dofs_l, V)
bc_r = dirichletbc(ScalarType(vr), b_dofs_r, V)
bc_t = dirichletbc(ScalarType(vt), b_dofs_t, V)
bc_b = dirichletbc(ScalarType(vb), b_dofs_b, V)

bcs = [bc_l, bc_r, bc_t, bc_b]

# Variational problem
f_ = Constant(domain, ScalarType(f))
dt_ = Constant(domain, ScalarType(dt))

#D_ = Expression('D', D=D)

t = 0

# solve problem
xdmf = XDMFFile(domain.comm, "output.xdmf", "w")

# Step 1
D_ = Constant(domain, ScalarType(D))

F = (u - u_n)*v*dx + dt_*D_*dot(grad(u),grad(v))*dx + dt_*f_*v*dx
a = form(lhs(F))
L = form(rhs(F))

# Solve problem
u_ = Function(V)

# assemble bilinear matrix
A = assemble_matrix(a, bcs=bcs)
A.assemble()

# assemble linear vector
b = assemble_vector(L)

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

# write value
xdmf.write_mesh(domain)
xdmf.write_function(u_n, t)

t += dt

# re-define assembled vector
with b.localForm() as loc_b:
    loc_b.set(0)
assemble_vector(b, L)

apply_lifting(b, [a], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)

solver.solve(b, u_.vector)
u_.x.scatter_forward()
    
# update u_n
u_n.x.array[:] = u_.x.array[:]

xdmf.write_function(u_n, t)

# step 2
t += dt

# change diffusion constant
D_ = Constant(domain, ScalarType(D*0.1))

F = (u - u_n)*v*dx + dt_*D_*dot(grad(u),grad(v))*dx + dt_*f_*v*dx
a = form(lhs(F))
L = form(rhs(F))

A = assemble_matrix(a, bcs=bcs)
A.assemble()

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

# re-define assembled vector
with b.localForm() as loc_b:
    loc_b.set(0)
assemble_vector(b, L)

apply_lifting(b, [a], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)

# solve problem
solver.solve(b, u_.vector)
u_.x.scatter_forward()
    
# update u_n
u_n.x.array[:] = u_.x.array[:]

xdmf.write_function(u_n, t)

xdmf.close()

dokken · April 25, 2022, 5:00am

See for instance: Hyperelasticity — FEniCSx tutorial
where T (Only the z-value) is updated within the time loop without redefining the form by calling:

    T.value[2] = n * tval0

I.e:
D_.value[0]= 0.1*D

igarashi.takahiro · April 25, 2022, 5:53am

Thank you for your reply.
I changed the part you pointed out and calculated again, but the results were different from that of the program where we re-defined the form shown above.
What is wrong with the program below?

# Module
from mpi4py import MPI
from dolfinx.mesh import create_rectangle
from dolfinx.fem import FunctionSpace, Function, Constant, locate_dofs_geometrical, dirichletbc, form, Expression
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import XDMFFile
from ufl import TrialFunction, TestFunction, dot, grad, dx, lhs, rhs
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import LinearProblem
import numpy as np

# Parameters
lx,ly,nx,ny = 10.0,5.0,10,5
vl,vr,vt,vb = 5.0,10.0,10.0,0.0

D = 1.0
f = 0.0
dt = 1.0
u0 = 0.0

# Mesh
domain = create_rectangle(MPI.COMM_WORLD, ([0, 0],[lx, ly]), [nx, ny])

# Function space
V = FunctionSpace(domain, ("CG", 1))

# Functions
u,v = TrialFunction(V),TestFunction(V)
u_n = Function(V)

# initial value
u_n.vector.array[:] = u0 * np.ones(u_n.vector.array.shape)

# Boundary condition
def left_boundary(x):
    return(np.isclose(x[0], 0.0))

def right_boundary(x):
    return(np.isclose(x[0], lx))

def top_boundary(x):
    return(np.isclose(x[1], ly))

def bottom_boundary(x):
    return(np.isclose(x[1], 0.0))

b_dofs_l = locate_dofs_geometrical(V, left_boundary)
b_dofs_r = locate_dofs_geometrical(V, right_boundary)
b_dofs_t = locate_dofs_geometrical(V, top_boundary)
b_dofs_b = locate_dofs_geometrical(V, bottom_boundary)

bc_l = dirichletbc(ScalarType(vl), b_dofs_l, V)
bc_r = dirichletbc(ScalarType(vr), b_dofs_r, V)
bc_t = dirichletbc(ScalarType(vt), b_dofs_t, V)
bc_b = dirichletbc(ScalarType(vb), b_dofs_b, V)

bcs = [bc_l, bc_r, bc_t, bc_b]

# Variational problem
f_ = Constant(domain, ScalarType(f))
dt_ = Constant(domain, ScalarType(dt))

#D_ = Expression('D', D=D)

t = 0

# solve problem
xdmf = XDMFFile(domain.comm, "output.xdmf", "w")

# Step 1
D_ = Constant(domain, ScalarType(D))

F = (u - u_n)*v*dx + dt_*D_*dot(grad(u),grad(v))*dx + dt_*f_*v*dx
a = form(lhs(F))
L = form(rhs(F))

# Solve problem
u_ = Function(V)

# assemble bilinear matrix
A = assemble_matrix(a, bcs=bcs)
A.assemble()

# assemble linear vector
b = assemble_vector(L)

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

# write value
xdmf.write_mesh(domain)
xdmf.write_function(u_n, t)

t += dt

# re-define assembled vector
with b.localForm() as loc_b:
    loc_b.set(0)
assemble_vector(b, L)

apply_lifting(b, [a], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)

solver.solve(b, u_.vector)
u_.x.scatter_forward()
    
# update u_n
u_n.x.array[:] = u_.x.array[:]

xdmf.write_function(u_n, t)

# step 2
t += dt

# change diffusion constant
D_.value = 0.1*D

# re-define assembled vector
with b.localForm() as loc_b:
    loc_b.set(0)
assemble_vector(b, L)

apply_lifting(b, [a], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)

# solve problem
solver.solve(b, u_.vector)
u_.x.scatter_forward()
    
# update u_n
u_n.x.array[:] = u_.x.array[:]

xdmf.write_function(u_n, t)

xdmf.close()

program (Re-defined form)

programReDefined1335×699 38.7 KB
program (New)

programNew1336×700 27.4 KB

dokken · April 25, 2022, 6:10am

Please start by simplifying your programs by using a for loop to avoid repeating code.

igarashi.takahiro · April 25, 2022, 6:55am

Sorry for the confusion.
I have rewritten it as follows using for loop.
We obtained the results with time dependent diffusion coefficient.
Thank you very much!

However, new question has emerged.
We will post that question to new thread.

# Module
from mpi4py import MPI
from dolfinx.mesh import create_rectangle
from dolfinx.fem import FunctionSpace, Function, Constant, locate_dofs_geometrical, dirichletbc, form, Expression
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import XDMFFile
from ufl import TrialFunction, TestFunction, dot, grad, dx, lhs, rhs
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import LinearProblem
import numpy as np

# Parameters
lx,ly,nx,ny = 10.0,5.0,10,5
vl,vr,vt,vb = 5.0,10.0,10.0,0.0

D = 1.0
f = 0.0
dt = 1.0
u0 = 0.0

# Mesh
domain = create_rectangle(MPI.COMM_WORLD, ([0, 0],[lx, ly]), [nx, ny])

# Function space
V = FunctionSpace(domain, ("CG", 1))

# Functions
u,v = TrialFunction(V),TestFunction(V)
u_n = Function(V)

# initial value
u_n.vector.array[:] = u0 * np.ones(u_n.vector.array.shape)

# Boundary condition
def left_boundary(x):
    return(np.isclose(x[0], 0.0))

def right_boundary(x):
    return(np.isclose(x[0], lx))

def top_boundary(x):
    return(np.isclose(x[1], ly))

def bottom_boundary(x):
    return(np.isclose(x[1], 0.0))

b_dofs_l = locate_dofs_geometrical(V, left_boundary)
b_dofs_r = locate_dofs_geometrical(V, right_boundary)
b_dofs_t = locate_dofs_geometrical(V, top_boundary)
b_dofs_b = locate_dofs_geometrical(V, bottom_boundary)

bc_l = dirichletbc(ScalarType(vl), b_dofs_l, V)
bc_r = dirichletbc(ScalarType(vr), b_dofs_r, V)
bc_t = dirichletbc(ScalarType(vt), b_dofs_t, V)
bc_b = dirichletbc(ScalarType(vb), b_dofs_b, V)

bcs = [bc_l, bc_r, bc_t, bc_b]

# Variational problem
f_ = Constant(domain, ScalarType(f))
dt_ = Constant(domain, ScalarType(dt))

#D_ = Expression('D', D=D)

t = 0

# solve problem
xdmf = XDMFFile(domain.comm, "output2.xdmf", "w")

# Step 1
D_ = Constant(domain, ScalarType(D))

F = (u - u_n)*v*dx + dt_*D_*dot(grad(u),grad(v))*dx + dt_*f_*v*dx
a = form(lhs(F))
L = form(rhs(F))

# Solve problem
u_ = Function(V)

# assemble bilinear matrix
A = assemble_matrix(a, bcs=bcs)
A.assemble()

# assemble linear vector
b = assemble_vector(L)

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

# write value
xdmf.write_mesh(domain)
xdmf.write_function(u_n, t)

for n in range(10):

    D_.value = 0.1*D*(n+1)

    t += dt

    # re-define assembled vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    assemble_vector(b, L)

    apply_lifting(b, [a], [bcs])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b, bcs)

    solver.solve(b, u_.vector)
    u_.x.scatter_forward()
    
    # update u_n
    u_n.x.array[:] = u_.x.array[:]

    xdmf.write_function(u_n, t)

xdmf.close()

Topic		Replies	Views
Speeding up time-dependent diffusion equation dolfinx	5	820	July 31, 2023
Diffusion behavior differs by changing constant value with .value method dolfinx	4	582	April 26, 2022
Dolfinx: Time-dependent problem equilibrating instantly dolfinx	1	550	April 6, 2022
Diffusion of a Gaussian function - source term on for fraction of simulation time dolfinx time-dependent	7	181	June 30, 2024
Time-dependent mesh and nonlinear problem (PETSC Error after migrating from FEniCSX 0.7.2 to 0.8.0) General	2	37	October 20, 2024

Time dependent diffusion constant in dolfinx

Related Topics