Diffusion behavior differs by changing constant value with .value method

dolfinx
1

Hello everyone.

The diffusion behavior differs between defining constant value with Constant method and changing constant value with .value method.
Below is a programs two patterns above.
Why does that happen?

Thank you for your attention.

  • define constant using donfinx.fem.Constant
# Module
from mpi4py import MPI
from dolfinx.mesh import create_unit_square
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
nx,ny = 10,10
vl,vr = 0.0,1.0

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

# Mesh
domain = create_unit_square(MPI.COMM_WORLD, 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], 1.0))

b_dofs_l = locate_dofs_geometrical(V, left_boundary)
b_dofs_r = locate_dofs_geometrical(V, right_boundary)

bc_l = dirichletbc(ScalarType(vl), b_dofs_l, V)
bc_r = dirichletbc(ScalarType(vr), b_dofs_r, V)

bcs = [bc_l, bc_r]

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

t = 0

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

D_ = Constant(domain, ScalarType(D*0.5))

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

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

画像_2022-04-25_164500

  • change constant using .value method
# Module
from mpi4py import MPI
from dolfinx.mesh import create_unit_square
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
nx,ny = 10,10
vl,vr = 0.0,1.0

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

# Mesh
domain = create_unit_square(MPI.COMM_WORLD, 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], 1.0))

b_dofs_l = locate_dofs_geometrical(V, left_boundary)
b_dofs_r = locate_dofs_geometrical(V, right_boundary)

bc_l = dirichletbc(ScalarType(vl), b_dofs_l, V)
bc_r = dirichletbc(ScalarType(vr), b_dofs_r, V)

bcs = [bc_l, bc_r]

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

t = 0

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

D_ = Constant(domain, ScalarType(D*1.0))

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

    t += dt

    D_.value = 0.5*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)

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

画像_2022-04-25_164532

2

You need to reassemble the matrix A as well, as it depends on D_

1 Like
3

Thank you for reply.
I tried to add the command for assembling matrix as you mentioned.
However, behavior was not changed with before.
Do I have any mistake?
Following is the program.

# Module
from mpi4py import MPI
from dolfinx.mesh import create_unit_square
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
nx,ny = 10,10
vl,vr = 0.0,1.0

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

# Mesh
domain = create_unit_square(MPI.COMM_WORLD, 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], 1.0))

b_dofs_l = locate_dofs_geometrical(V, left_boundary)
b_dofs_r = locate_dofs_geometrical(V, right_boundary)

bc_l = dirichletbc(ScalarType(vl), b_dofs_l, V)
bc_r = dirichletbc(ScalarType(vr), b_dofs_r, V)

bcs = [bc_l, bc_r]

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

t = 0

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

D_ = Constant(domain, ScalarType(D*1.0))

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.5*D

    # re-define assembled matrix
    A = assemble_matrix(a, bcs=bcs)
    A.assemble()

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

You shouldn’t assemble a new matrix, but rather reassemble into your existing A which has been assigned to the KSP with the solver.setOperators(A) call.

5

Thank you for reply.
After adding “solver.setOperators(A)” after “A.assemble()”, I got the expected result. Thank you very much!
Do I have any mistake about reassembling operation of matrix in following program?

# Module
from mpi4py import MPI
from dolfinx.mesh import create_unit_square
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
nx,ny = 10,10
vl,vr = 0.0,1.0

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

# Mesh
domain = create_unit_square(MPI.COMM_WORLD, 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], 1.0))

b_dofs_l = locate_dofs_geometrical(V, left_boundary)
b_dofs_r = locate_dofs_geometrical(V, right_boundary)

bc_l = dirichletbc(ScalarType(vl), b_dofs_l, V)
bc_r = dirichletbc(ScalarType(vr), b_dofs_r, V)

bcs = [bc_l, bc_r]

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

t = 0

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

D_ = Constant(domain, ScalarType(D*1.0))

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.5*D
    A = assemble_matrix(a, bcs=bcs)
    A.assemble()
    solver.setOperators(A)

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