Error when trying to implement continuation method to run fluid flow problem at higher Reynolds numbers

General
,
1

Hi,

I am trying to implement the method called ‘continuation’ discussed below, using a previous solution at some viscosity, to get a solution for a lower viscosity, by using the higher viscosity solution as an initial guess for the solution.

In the previous version of Fenics this was done using the assign() function, now it gives an error. What could I use now instead of that function?

I also tried:

u.vector()[:] = u_old.vector()
p.vector()[:] = p_old.vector()

...

u_old.vector()[:] = u.vector()
p_old.vector()[:] = p.vector()

A minimal running code giving the error is below.

Thank you,

Alex

from dolfinx import mesh, fem, io, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

from mpi4py import MPI

import numpy as np

import ufl
from ufl import inner, div, grad, dx

import basix

import time

from petsc4py import PETSc


height = 0.06
width  = 0.08

n_h = 20
n_w = 20


msh = mesh.create_box(MPI.COMM_WORLD, ((-width/2, -width/2, 0.0), (width/2,width/2,height)), (n_h,n_h,n_w), mesh.CellType.tetrahedron)

fdim = msh.topology.dim - 1

rho_air = fem.Constant(msh, 1.225)


# Function to mark noslip
def noslip_boundary1(x):
    return np.isclose(x[2], height)


P2 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 2, shape=(msh.geometry.dim, ))
P1 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 1)

V, Q = fem.functionspace(msh, P2), fem.functionspace(msh, P1)

# Create the function space
TH = basix.ufl.mixed_element([P2, P1])
W = fem.functionspace(msh, TH)
W0, _ = W.sub(0).collapse()


# No slip boundary condition 1
noslip  = fem.Function(V)
facets1 = mesh.locate_entities_boundary(msh, fdim, noslip_boundary1)
dofs1   = fem.locate_dofs_topological((W.sub(0), V), fdim, facets1)
bc0     = fem.dirichletbc(noslip, dofs1, W.sub(0))


# Since for this problem the pressure is only determined up to a constant, we pin it at a point
zero = fem.Function(Q)
zero.x.array[:] = 0
dofs = fem.locate_dofs_geometrical(
    (W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))



def body_force_fun(x):
    F_max    = 2.0e-2
    z_offset = 0.02

    r_xy = np.sqrt(x[0]**2 + x[1]**2)

    a = 0.005
    max_val = a * np.exp(-0.5)
    Fr_val  = np.exp(-r_xy**2/(2*a**2)) / max_val
        
    b = 0.005
    Fz_val = np.exp(-(x[2]-z_offset)**2/(2*b**2))
    F_mag  = F_max * Fr_val * Fz_val
      
    # -ve sign to get clockwise rotation from array top
    return -F_mag * np.vstack((x[1], -x[0], np.zeros_like(x[2])))
        

Force_per_unit_vol = fem.Function(V)
Force_per_unit_vol.interpolate(lambda x: body_force_fun(x))




# Collect Dirichlet boundary conditions
bcs = [bc0, bc2]

# Define variational problem
w     = fem.Function(W)
w_old = fem.Function(W)

(u, p) = ufl.split(w)
(u_old, p_old) = ufl.split(w_old)

(v, q) = ufl.TestFunctions(W)



for mair in [4.0e-5, 1.8e-5]:

    mu_air  = fem.Constant(msh, mair)

    u.assign(u_old)
    p.assign(p_old)

    F = rho_air*inner(grad(u)*u, v)*dx + \
        mu_air*inner(grad(u), grad(v))*dx - \
        inner(Force_per_unit_vol, v)*dx - \
        inner(p, div(v))*dx - \
        div(u)*q*dx
    
    dw = ufl.TrialFunction(W)
    dF = ufl.derivative(F, w, dw)
    
    problem = NonlinearProblem(F, w, bcs=bcs, J=dF)
    
    solver = NewtonSolver(MPI.COMM_WORLD, problem)
    solver.convergence_criterion = "incremental"
    solver.rtol   = 1e-6
    solver.report = True
    solver.error_on_nonconvergence = False
    
    
    # Custom solver options 
    opts = PETSc.Options()  # type: ignore
    ksp = solver.krylov_solver
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}ksp_type"] = "preonly"
    opts[f"{option_prefix}pc_type"] = "lu"
    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
    ksp.setFromOptions()
    
    
    # Compute the solution
    log.set_log_level(log.LogLevel.INFO)
    start = time.time()
    solver.solve(w)
    end   = time.time()
    
    u_old.assign(u)
    p_old.assign(p)
    
    print('\n')
    print(msh.topology.index_map(msh.topology.dim).size_local)
    print((end - start) / 60)
    print('\n\n')



# Split the mixed solution and collapse
u = w.sub(0).collapse()
p = w.sub(1).collapse()

# Write the solution to file
with io.XDMFFile(MPI.COMM_WORLD, "pressure.xdmf", "w") as pfile_xdmf:
    p.x.scatter_forward()
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(p)

with io.XDMFFile(MPI.COMM_WORLD, "velocity.xdmf", "w") as ufile_xdmf:
    u.x.scatter_forward()
    P1 = basix.ufl.element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
    u1 = fem.Function(fem.functionspace(msh, P1))
    u1.interpolate(u)
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u1)
    
with io.XDMFFile(MPI.COMM_WORLD, "body_force.xdmf", "w") as ufile_xdmf:
    Force_per_unit_vol.x.scatter_forward()
    P1 = basix.ufl.element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
    u1 = fem.Function(fem.functionspace(msh, P1))
    u1.interpolate(Force_per_unit_vol)
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u1)
2

fe.assign and u.vector()[:] has been replaced with function.x.array[:].

So in your loop, replace the fe.assign for (u, p) with

w.x.array[:] = w_old.x.array[:]

Similarly, at the end

w_old.x.array[:] = w.x.array[:]
3

Thank you Nick, this works