Assigning value to a Function within the time stepping loop

dolfinx
1

Hi,
I am trying to implement the legacy FEinCS’s advection, diffusion, reaction code using dolfinx. In particular, I am trying to assign values to a vector function w (velocity - previously generated by the Navier Stokes code) by reading data from a stored xdmf file. As this w is used in the variational formulation, I suppose I need to reassemble the form, F.

I am aware of the checkpointing issue, therefore, have tried meshio in the code below.

Here is my failed attempt!
Any help is deeply appreciated!

Thanks in advance.

"""
Trying the generalised reaction, diffusion in fenicsx
"""
import numpy as np
import tqdm.autonotebook
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.fem import (Constant, Function, FunctionSpace, VectorFunctionSpace, petsc,
                         assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                               create_vector, create_matrix, set_bc, NonlinearProblem)
from dolfinx.nls.petsc import NewtonSolver
from ufl import (FacetNormal, FiniteElement, MixedElement, Identity, Measure,
                 TestFunctions, TrialFunctions, TrialFunction, TestFunction,
                 VectorElement, split, SpatialCoordinate,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
from dolfinx.io import XDMFFile
from dolfinx import plot
from dolfinx import nls
from dolfinx.io.gmshio import read_from_msh
import pyvista as pv
import meshio
from petsc4py.PETSc import ScalarType

# Load external mesh domain that was generated using GMSH
domain, cell_tags, ft = read_from_msh(filename="mesh2d.msh",
                                      comm=MPI.COMM_WORLD,
                                      rank=0,
                                      gdim=2)
gdim = domain.topology.dim
fdim = gdim - 1

# Program
t = 0
T = 3.0
num_steps = 1500
dt = T / num_steps
gamma = 0.1
K = 0.1

# Function space
P1 = FiniteElement(family="Lagrange", cell=domain.ufl_cell(), degree=1)
V = FunctionSpace(mesh=domain, element=P1)

v1 = TestFunction(function_space=V)
u1 = Function(V=V)
u1n = Function(V=V)


# Creating Functionspace for velocity data
# Pvec = VectorElement(family="Lagrange", cell=domain.ufl_cell(), degree=2)
Pvec = VectorElement(family="CG", cell=domain.ufl_cell(), degree=2)
W = FunctionSpace(mesh=domain, element=Pvec)

w = Function(V=W)


# Source function
class F1_Expression:
    def __init__(self):
        self.t = 0.0

    def eval(self, x):
        values = np.full(x.shape[1], 0.0)
        idx = (x[0] - 0.25) ** 2 + (x[1] - 0.75) ** 2 < 0.05 ** 2
        values[idx] = 0.1
        return values


# inheriting from class definition
f1_expression = F1_Expression()
f1_expression.t = 0

# interpolating into a function
f1 = Function(V=V)
f1.interpolate(u=f1_expression.eval)

# defining the constants to be used in the weak form
kdt = Constant(domain=domain, c=PETSc.ScalarType(dt))
K = Constant(domain=domain, c=PETSc.ScalarType(K))
gamma = Constant(domain=domain, c=PETSc.ScalarType(gamma))

# Define variational problem
F = ((u1 - u1n) / kdt) * v1 * dx + dot(w, grad(u1)) * v1 * dx +\
    gamma * dot(grad(u1), grad(v1)) * dx + K * u1 * v1 * dx - f1 * v1 * dx


# Problem solver configuration
problem = NonlinearProblem(F, u1)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
# solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True


# Loading the velocity data
reader = meshio.xdmf.TimeSeriesReader("w.xdmf")
reader.read_points_cells()  # Necessary


# Class to set value of w
class SetValueOfW:
    def __init__(self, t, data):
        self.t = t
        self.data = data

    def eval(self, x):
        return np.ones(self.data.shape)


w_value = SetValueOfW(t=t, data=np.transpose(reader.read_data(k=0)[1]['u'][:, 0:2]))
print("t={}".format(t))


# Time stepping and solving
for n in range(num_steps):
    t += dt
    # w.x.array = reader.read_data(k=n)[1]['u']
    print("t={}".format(t))

    # I suppose I am correctly reading and assigning the values to w
    # but the change fails to reflect in the time loop.
    w_value.t = t
    w_value.data = np.transpose(reader.read_data(k=n)[1]['u'][:, 0:2])
    w.interpolate(w_value.eval)

    r = solver.solve(u1)
    print(f"Step {int(t / dt)}: num iterations: {r[0]}")
    u1n.x.array[:] = u1.x.array


print("Done!")

Here is the message from the solver.

  File "/Users/anil/Apps/anaconda3/envs/fenics-env/lib/python3.8/site-packages/dolfinx/nls/petsc.py", line 41, in solve
    n, converged = super().solve(u.vector)
RuntimeError: Newton solver did not converge because maximum number of iterations reached