PETSc error while working with mixed form


I am solving time dependent fluid flow poblem in dolfin (2019.01) and I am struggling between interpolation and PETSc error while working with mixed elements. If I manage to get rid of the interpolation error then PETSc is unhappy while fixing the PETSc causes the interpolation become impossible. Additionally, I need to solve it using nonlinear solver instead of the linear algorithms illustrated in tutorials. I will sincerely appreciate help from the experts.

Below is an MWE:

import gmsh, meshio, dolfin
import matplotlib.pyplot as plt
import tqdm.autonotebook
import numpy as np
from scipy.optimize import fsolve

# Print log messages only from the root process in parallel
dolfin.parameters["std_out_all_processes"] = False;
# Enable extrapolation when mapping fields between different meshes
dolfin.parameters['allow_extrapolation'] = True

def generate_gmsh_mesh():
    gdim = 2
    model_rank = 0
    occ = gmsh.model.occ
    fov = (1, 1)
    rectangle = occ.addRectangle(0,0,0, fov[0], fov[1])

    domains = gmsh.model.getEntities(gdim)
    for j, dom in enumerate(domains):
        gmsh.model.addPhysicalGroup(gdim, [dom[1]], tag=j+1)

    all_bnd = gmsh.model.getEntities(gdim - 1)
    for bnd in all_bnd:
        gmsh.model.addPhysicalGroup(gdim - 1, [bnd[1]], tag=bnd[1])


    # close gmsh

    # obtain dolfin mesh
    msh ="dolfin_ss_mesh.msh")
    ############### new block for saving and loading mesh with markers
    def create_mesh(mesh, cell_type, prune_z=True):
        cells = mesh.get_cells_type(cell_type)
        cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
        mesh_points = mesh.points[:, :2] if prune_z else mesh.points
        out_mesh = meshio.Mesh(points=mesh_points, 
                               cells={cell_type: cells}, 
                               cell_data={"name_to_read": [cell_data]})
        return out_mesh
    triangle_mesh = create_mesh(msh, "triangle",prune_z=True)
    line_mesh = create_mesh(msh, "line",prune_z=True)
    meshio.write("dolfin_ss_mesh.xdmf", triangle_mesh)
    meshio.write("dolfin_ss_mesh_mf.xdmf", line_mesh)

    # now read it from xdmf format
    mesh = dolfin.Mesh()
    with dolfin.XDMFFile("dolfin_ss_mesh.xdmf") as infile:
    mvc = dolfin.MeshValueCollection("size_t", mesh, dim=1) 
    with dolfin.XDMFFile("dolfin_ss_mesh_mf.xdmf") as infile:, "name_to_read")
    mesh_bnd = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    return mesh, mesh_bnd

dt, T = 0.01, 1
num_steps = dt/T
Um = 1
L = 1
rhof, mu = 1, 1

u_in = dolfin.Expression(("4*{}*x[1]*({}-x[1])/{}".format(Um, L, L**2), "0.0"), degree=2)
u_noslip = dolfin.Constant((0.0, 0.0))
p_out = dolfin.Constant(0.0)

# initialize mesh
mesh, mesh_bnd = generate_gmsh_mesh()

# Define coefficients
k = dolfin.Constant(dt)
f = dolfin.Constant((0, 0))

inner = dolfin.inner
dx = dolfin.dx
grad = dolfin.grad
div = dolfin.div

V_el = dolfin.VectorElement("CG", mesh.ufl_cell(), 2)
Q_el = dolfin.FiniteElement("CG", mesh.ufl_cell(), 1)
W = dolfin.FunctionSpace(mesh, V_el*Q_el)
up0 = dolfin.Function(W)
u0 = up0.split()[0] # get a Function

# Time-stepping
progress = tqdm.autonotebook.tqdm(desc="Solving PDE", total=num_steps)
t = dt # time for the first iteration
t_idx = 1
while t < T + dolfin.DOLFIN_EPS:
    # generate new mesh given the position of the cylinder
    mesh, mesh_bnd = generate_gmsh_mesh() # in MWE, the mesh does not change, but later it will
    # define the physics given the new mesh

    # Define function spaces (P2-P1)
    V_el = dolfin.VectorElement("CG", mesh.ufl_cell(), 2)
    Q_el = dolfin.FiniteElement("CG", mesh.ufl_cell(), 1)
    W = dolfin.FunctionSpace(mesh, V_el*Q_el)
    ut, pt = dolfin.TrialFunctions(W)
    up = dolfin.Function(W)
    upn = dolfin.Function(W)
    u, p = up.split(deepcopy=True)
    un, pn = upn.split(deepcopy=True) 
    u.interpolate(u0) # set the starting point of the nonlinear solver to the old solution
    un.interpolate(u0) # interpolate old solution to the new mesh
    F = rhof/k*inner(u - un, ut)*dx + rhof*inner(grad(u)*u, ut)*dx + mu*inner(grad(u), grad(ut))*dx - p*div(ut)*dx 
    + rhof*pt*div(u)*dx

    x = dolfin.SpatialCoordinate(mesh)
    # boundary conditions   

    bc_bot = dolfin.DirichletBC(W.sub(0), u_noslip, mesh_bnd, 1)
    bc_outlet = dolfin.DirichletBC(W.sub(1), p_out, mesh_bnd, 2)
    bc_top = dolfin.DirichletBC(W.sub(0), u_noslip, mesh_bnd, 3)
    bc_inlet = dolfin.DirichletBC(W.sub(0), u_in, mesh_bnd, 4)

    bcs = [bc_top, bc_bot, bc_inlet, bc_outlet]
    dolfin.solve(F==0, up, bcs=bcs, 
                 solver_parameters={"newton_solver":{"relative_tolerance": 1e-6,
                                                     # "linear_solver": "mumps", 
                                                     "maximum_iterations" : 100,
                 form_compiler_parameters={"cpp_optimize": True, "representation": "uflacs"} )
    # save the current distribution of velocity and pressure
    ufile << u
    pfile << p
#     Prepare for the next time step
    up0 = dolfin.Function(W) # Function for the current mesh
    u0 = up0.split(deepcopy=True)[0]

    u0.assign(u) # store the current solution for the next iteration
    t += dt
    t_idx += 1

I suspect my issue is related to dolfin.split(u) vs u.split()

After searching through threads on this forum, I am now firm that for defining the variational form, I need
u, p = dolfin.split(up) # (right?)

As for interpolation, I am not sure the right way since dolfin.split(u) returns a ListTensor object which does not have an interpolation method. I came across this though:

Following that, here is how I should do it in the loop (the rest remains unchanged):

V_el = dolfin.VectorElement("CG", mesh.ufl_cell(), 2)
Q_el = dolfin.FiniteElement("CG", mesh.ufl_cell(), 1)
W = dolfin.FunctionSpace(mesh, V_el*Q_el)
up = dolfin.Function(W) # current step solution
upn = dolfin.Function(W) # previous step solution 
u, p = dolfin.split(up) # for defining the variational form
un, pn = dolfin.split(upn) # for defining the variational form

up.sub(0).interpolate(u0) # set initial starting point for newton solver

However, this now leads me to the following error:

Any help will be highly appreciated

To set initial conditions to a Mixed space you have to use the FunctionAssigner, Ref Usage of FunctionAssigner - #2 by dokken

Sorry I do not get how the example translates to my case. In the referred link, two functions defined on a scalar space are being assigned to a function defined on a vector space.

In my case, both up and up0 are functions out of the same mixed element space. u and u0 are the first split components of up and up0. Could you kindly show how the assignment step will be written in my case?

Thanks a lot!

There is no significant different between assigning scalar components from a scalar space to a vector space and a space to a vector space, as shown in:

The documentation is straightforward: and there are examples: pyadjoint/ at e1c2292a92398f9a558b03495dc8a774a2d94bd0 · dolfin-adjoint/pyadjoint · GitHub