How to store mixed solution from Function and set it as initial guess for another solution

General
1

Hi,

This question is related to one I posted previously: Error when trying to implement continuation method to run fluid flow problem at higher Reynolds numbers

That has running code. In the legacy Fenics this could be done with:

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

Based on some of the answers here, it seems one way to do it is using:

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

But I get the following error:

'ListTensor' object has no attribute 'x'

How is the assign() function replaced in FEniCSx?

Thank you,

Alex

2

I suppose you used u_old,p_old = ufl.split(solution) to obtain u_old and p_old?

It depends a little on how you set it all up (hence a MWE is important), but you can probably do what you want to do by:

u_old,p_old = [ subsol.collapse() for subsol in solution.split() ]

Note this not only splits the solution into it subcomponents, it also collapses them into their stand-alone representations.

If you need something that avoids recreating the u_old and p_old function objects then you’d have to supply a MWE.

3

Hi Stein,

Thank you for the prompt reply. A MWE is shown below. In the code I’m trying to use a previous solution to solve the same fluid flow problem at a lower viscosity.

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, "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)
4

As you’re always assigning the complete (mixed) solution, wouldn’t you just be able to:

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

Thank you that worked!

It seems to help with convergence but not by a lot. Do you know any other techniques (besides using a turbulence model) to be able to run Fenics Navier-Stokes code for higher Reynolds numbers?

6

Almost by definition of the term “turbulence model”, I am afraid the answer is no.

As the Reynolds number increases the PDE becomes more unstable. Past its critical value turbulence kicks in, characterized by a wide range in spatiotemporal scales. On a computational grid of managable size (i.e., not DNS resolved), this simply means you cannot represent the true solution to the PDE. But maybe you don’t have to; a coarse approximation to the true solution (=LES) is typically all you need to obtain your quantities of interest.

But, recall that our FEM error is a combination of interpolation error (\approx difference between DNS and optimal LES) and discretization/Galerkin error. So, in order to get a good LES solution, the latter needs to be small. That’s where you have to deal with the unstable nature of the PDE through stabilization, and that’s honestly all turbulence models do.

In FEM, we typically stabilize advection based PDE’s with SUPG. It turns out, from the Variational Multiscale perspective, that such a stabilized method can be derived as a scale-interaction model = a turbulence model (see Redirecting).

So, TL;DR: in FEM you can use ‘regular’ stabilization methods as LES turbulence models.

7

That’s a great answer thank you! I did take class on turbulence in my masters and don’t remember learning about this. It actually makes more sense. Is this one example of what you describe?

Do you know of any Fenics code that uses a simpler stabilization method?

8

I suppose they do. Their equations on the bottom of section 3.1 are very similar to those in the work I cited ealier. I’m not sure why they don’t use the usual nomenclature for describing that as SUPG PSPG and LSIC (sometimes “grad-div”) stabilization.

On the forum you’ll find various examples of people having implemented SUPG for Navier-Stokes. In combination with Taylor-Hood elements, that should get you pretty far.