Multiple equation system by individually iterating on each equation

1

Hello,

I’m trying to solve a multiple equation system by using the Newton solver on each equation in the system iteratively, similar to how pressure, velocity and turbulence model corrections are handled in separate correction steps for the Reynolds-averaged Navier-Stokes equations in solvers like simpleFoam in OpenFOAM.

I have been trying to set up a simple example of a set of equations solved in this manner by adapting the convection-diffusion equation tutorial, but without much success so far. Specifically, I have been trying to set up a system of two coupled, steady, nonlinear convection-diffusion equations in 2D and solve them by

  • Create initial guesses for both concentrations c1 and c2
  • In a loop, iterate until convergence:
    1. Plug the current guess for c2 into equation 1 and solve for c1
    2. Plug the current guess for c1 into equation 2 and solve for c2

This is the current code I use for this purpose:

import fenics
import numpy as np
    class SteadyConvectionDiffusionGenerator:
        def __init__(self):
            self.ndims = 2

        #mesh
        p_bottomleft = fenics.Point(*[-2.0 for _ in range(self.ndims)])
        p_topright = fenics.Point(*[2.0 for _ in range(self.ndims)])
        self.mesh = fenics.RectangleMesh.create(p=[p_bottomleft,p_topright],n=[100 for _ in range(self.ndims)],cell_type=fenics.CellType.Type.triangle)

        #finite element setup stuff
        P1 = fenics.FiniteElement('CG', fenics.triangle, 1)
        element = P1
        self.FuncSpace = fenics.FunctionSpace(self.mesh, element)
        self.TestFuncs = [fenics.TestFunctions(self.FuncSpace) for _ in range(self.ndims)]
        self.SolnFuncs = [fenics.Function(self.FuncSpace) for _ in range(self.ndims)]
        self.u = [fenics.split(sf) for sf in self.SolnFuncs]

        #forcing
        f_1 = fenics.Expression('pow(2*x[0]-0.1,2)+pow(x[1]-0.1,2)<0.05*0.05 ? 0.1 : 0',degree=1)
        f_2 = fenics.Expression('pow(2*x[0]-0.1,2)+pow(x[1]-0.3,2)<0.05*0.05 ? 0.1 : 0',degree=1)
        
        #equations
        diffusion_coeff = fenics.Constant(0.01)
        closure_coeff = fenics.Constant(10.0)

        self.eq1_linear_part = \
            + fenics.inner(fenics.grad(self.u[0][0]),fenics.grad(self.TestFuncs[0][0]))*diffusion_coeff*fenics.dx \
            - f_1 * self.TestFuncs[0][0] * fenics.dx
        self.eq1 = self.eq1_linear_part + closure_coeff * self.u[0][0] * self.u[1][0] * self.TestFuncs[0][0] * fenics.dx
            
        self.eq2_linear_part = \
            + fenics.inner(fenics.grad(self.u[1][0]),fenics.grad(self.TestFuncs[1][0]))*diffusion_coeff*fenics.dx \
            - f_2 * self.TestFuncs[1][0] * fenics.dx
        self.eq2 = self.eq2_linear_part + closure_coeff * self.u[0][0] * self.u[1][0] * self.TestFuncs[1][0] * fenics.dx

        self.sp = {
            "newton_solver":{"maximum_iterations":50, 'relative_tolerance': 1e-6, 'absolute_tolerance': 1e-5, "error_on_nonconvergence": False, "linear_solver": "lu"}
        }

    def solve_together(self):
        fenics.solve(self.eq1 + self.eq2 == 0, self.SolnFuncs[0], solver_parameters = self.sp)

    def solve_one_by_one(self):
        for k in range(10):
            print('---Solving eq 1---')
            fenics.solve(self.eq1 == 0, self.SolnFuncs[0], solver_parameters = self.sp)
            print('---Solving eq 2---')
            fenics.solve(self.eq2 == 0, self.SolnFuncs[1], solver_parameters = self.sp)
            print('\n')

if __name__ == '__main__':
    
    g = SteadyConvectionDiffusionGenerator()
    # g.solve_together()
    g.solve_one_by_one()

This code runs, but I am unable to export or plot the solutions to see if the solutions are as expected, as I get the following error for the VTK export functionality:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to write data to VTK file.
*** Reason: Don’t know how to handle vector function with dimension other than 2 or 3.
*** Where: This error was encountered inside VTKFile.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

However, when I inspect the solution I get by looking at the raw values in g.SolnFuncs[0].vector() for example, I see that they are on the order of 1e+12 which should not be the case for this problem, so there’s certainly something wrong here.

I could not find much info on how to do this kind of thing with Fenics - there’s this unanswered post from about a year ago: How to decouple and iteratively solve Cahn-Hilliard equation

First of all, how can I fix the export problem? Secondly, could anyone point me towards resources showing how to implement such a solution algorithm?

2

For the first issue, you can use XDMFFile instead i of File (VTK/pvd) to save your solution. Note that you should use write_checkpoint if you want to use other spaces than CG 1, see: Loading xdmf data back in

This might help you to debug your problem

3

XDMF seems to have fixed the file export problem, thanks for the recommendation. Now, I see that the solution values are fairly wrong and are largely uniform over the solution domain - one variable has constant values of around 3e+12 and the other around 2e-27. Any clues as to why this may be?

4

I have not had time to look at your code. A thing that strikes me as Odd is the following

Why are you creating a Mixed element with only one element?

Does the coupled solver Give the same as the sequential one?

5

Thanks for pointing that out - I had actually fixed that little error but forgot to paste the updated version. I will rectify it in the original post as well. This should be fine in the sequential version I hope:

P1 = fenics.FiniteElement('CG', fenics.triangle, 1)
element = P1

The coupled version of the code converges when I use an approach similar to the one in the tutorial. The code for that is as follows:

class SteadyConvectionDiffusionGenerator:
    def __init__(self, ndims = 2):
        self.ndims = 2

        #mesh
        p_bottomleft = fenics.Point(*[-2.0 for _ in range(ndims)])
        p_topright = fenics.Point(*[2.0 for _ in range(ndims)])
        self.mesh = fenics.RectangleMesh.create(p=[p_bottomleft,p_topright],n=[100 for _ in range(ndims)],cell_type=fenics.CellType.Type.triangle)

        #finite element setup stuff
        P1 = fenics.FiniteElement('CG', fenics.triangle, 1)
        element = fenics.MixedElement([P1 for _ in range(ndims)])
        self.FuncSpace = fenics.FunctionSpace(self.mesh, element)
        self.TestFuncs = fenics.TestFunctions(self.FuncSpace)
        self.SolnFuncs = fenics.Function(self.FuncSpace)
        self.u = fenics.split(self.SolnFuncs)

        #forcing
        f_1 = fenics.Expression('pow(2*x[0]-0.1,2)+pow(x[1]-0.1,2)<0.05*0.05 ? 0.1 : 0',degree=1)
        f_2 = fenics.Expression('pow(2*x[0]-0.1,2)+pow(x[1]-0.3,2)<0.05*0.05 ? 0.1 : 0',degree=1)
        
        #equations
        diffusion_coeff = fenics.Constant(0.01)
        closure_coeff = fenics.Constant(10.0)

        self.eq1_linear_part = \
            + fenics.inner(fenics.grad(self.u[0]),fenics.grad(self.TestFuncs[0]))*diffusion_coeff*fenics.dx \
            - f_1 * self.TestFuncs[0] * fenics.dx
        self.eq1 = self.eq1_linear_part + closure_coeff * self.u[0] * self.u[1] * self.TestFuncs[0] * fenics.dx
            
        self.eq2_linear_part = \
            + fenics.inner(fenics.grad(self.u[1]),fenics.grad(self.TestFuncs[1]))*diffusion_coeff*fenics.dx \
            - f_2 * self.TestFuncs[1] * fenics.dx
        self.eq2 = self.eq2_linear_part + closure_coeff * self.u[0] * self.u[1] * self.TestFuncs[1] * fenics.dx
            

        #initialize soln randomly otherwise nan values are obtained
        self.SolnFuncs.vector().set_local(np.linspace(0.0,1.0,self.SolnFuncs.vector().size()))
        
        self.SolnFuncs.vector().apply("")

        self.sp = {
            "newton_solver":{"maximum_iterations":50, 'relative_tolerance': 1e-6, 'absolute_tolerance': 1e-5, "error_on_nonconvergence": False, "linear_solver": "lu"}
        }

    def solve_together(self):
        fenics.solve(self.eq1 + self.eq2 == 0, self.SolnFuncs, solver_parameters = self.sp)
            
if __name__ == '__main__':
    
    g = SteadyConvectionDiffusionGenerator()
    g.solve_together()
6

If I understood correctly, you are solving system of equations in a segregated manner. For iillustration of this kind of problem, please check this example. It is possible that the example is outdated, but you can probably figure it out how to adapt it for your problem.

7

Thanks, I’ll study that example in detail. Meanwhile, are you aware of the application of the method in this tutorial to cases where one or more sub-problems are nonlinear?

8

Solving of the system of nonlinear advection–diffusion–reaction equations is illustrated in this example. The problem is solved in coupled manner, but combining it with the previous example, it should be straightforward to implement it in segregated manner. I think that updated versions of both codes can be found here, in documented and undocumented demos.

9

I was aware of these tutorials and was following them, but I suppose I’ll spend a bit more time to make it work with a different approach so that each sub-problem is a linear one as the naive approach with the Newton solver seems to not converge.

That being said, are there any examples/tutorials on the usage of Fenics for the Reynolds-averaged Navier-Stokes (RANS) equations? Went through all the tutorials but could not find any examples of RANS.

10

Unfortunately, I do not know if there is tutorial or example for these equations.

Regarding your original problem, maybe you can try something like this:

V = FunctionSpace(mesh, 'P', 1)
u1 = TrialFunction(V)
v1 = TestFunction(V)
u2 = TrialFunction(V)
v2 = TestFunction(V)
u1_solution = Function(V)
u2_solution = Function(V)

F1 = eq1_linear_part + closure_coeff * u1 * u2_solution * v1 * dx # define variational form for the first equation, figuring solution of the second equation
F2 = eq2_linear_part + closure_coeff * u1_solution * u2 * v2 * dx # define variational form for the second equation, figuring solution of the first equation

F1 = action(F1, u1_solution)
F2 = action(F2, u2_solution)

u1_solution = ... # some initial guess
u2_solution = ... # some initial guess

for k in range(10):
	solve(F1 == 0, u1_solution,  bc, solver_parameters={"newton_solver": {"relative_tolerance": 1e-4}})
	solve(F2 == 0, u2_solution,  bc, solver_parameters={"newton_solver": {"relative_tolerance": 1e-4}})

I used something similar for the time dependent drift-diffusion equation, before I switch to coupled approach.