Multivariate Nonlinear Function

variational formulation
1

I have a system with two weak PDEs of the form:
F(u,v,q1)=0
and
G(u,v,q2)=0

I want to solve F for u with a temporary v, then vice versa for G.

I’m struggling with how to implement this.

The only way I can get things to run at the moment is to have u,v as Functions(V) and q1,q2 as q1,q2=TestFunctions(V).
But then it’s unclear what the solutions to NonVariationalProblem(F,x) are, is x an estimate for u or v? I thought it made the most sense to have a

u=TrialFunction(V)
v=TrialFunction(V)
tmpu=Function(V)
tmpv=Function(V)

but this gives me

File "/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py", line 41, in nonlinear_operator
    raise ArityMismatch("Applying nonlinear operator {0} to expression depending on form argument {1}.".format(o._ufl_class_.__name__, t))

My real question is, if the NonVariationalProblem contains two function terms then what does the solution represent?

MWE:



from dolfin import *
import random
import numpy as np


def Initial(perm,x,out):
    size=out.vector().size()
    r1=np.random.rand(size)
    r2=np.random.rand(size)
    out.vector().set_local(x+perm*(r1-r2))
    return out
def lapla(u):
    return inner(grad(u),grad(u))

def Simulations():

    # nx,ny are number of elements on x and y axis
    # xm, ym are maximum x and y values
    # x0,y0 are minimum x and y values
    # hx, hy are meshwidth on the x and y axis.
    x0=0
    y0=0
    xm=10
    ym=10
    nx=10
    ny=10
    StSt=[0,1]
    perm=0.1

    # FeNiCs variables--------------------------------------------------
    # Initialise the space ---------------------------------------------
    mesh=RectangleMesh(Point(x0,y0),Point(xm,ym),nx,ny)
    #P1 = FiniteElement('P', triangle, 1)
    #element = MixedElement([P1, P1])
    V = FunctionSpace(mesh, 'Lagrange',1)
    # Function(V) indicates a scalar variable over the domain.
    # TestFunction(V) is a scalar test function over the domain for the weak problem.
    u0    = TrialFunction(V)
    v0    = TrialFunction(V)
    tmpu0 = Function(V)
    tmpv0 = Function(V)
    u1    = Function(V)
    v1    = Function(V)
    q1      = TestFunction(V)
    q2      = TestFunction(V)



    #Initialise time step 0
    tmpu0 = Initial(perm,StSt[0],tmpu0)
    tmpv0 = Initial(perm,StSt[1],tmpv0)


    # Define variational problem ---------------------------------------
    F=(lapla(u0)+u0**2-u0*tmpv0)*q1*dx
    G=(lapla(v0)+(tmpu0+v0)*(1-v0)-v0)*q2*dx

    # Solvers
    problem_u=NonlinearVariationalProblem(F,u1)
    solver_u=NonlinearVariationalSolver(problem_v)
    tmpu0.assign(u1)
    problem_v=NonlinearVariationalProblem(G,v1)
    solver_v=NonlinearVariationalSolver(problem_v)
    tmpv0.assign(v1)
    return

def main(simon=1):
    if simon: Simulations()
    return

if __name__=='__main__':
    main()
    exit()
2

See, for example, the nonlinear Poisson demo. You’re casting the nonlinear problem with TrialFunctions which does not make sense.

There are actually a few things about your code which are very puzzling, in particular your variational formulations… despite this, the following code runs:

from dolfin import *
import random
import numpy as np


def Initial(perm,x,out):
    size=out.vector().size()
    r1=np.random.rand(size)
    r2=np.random.rand(size)
    out.vector().set_local(x+perm*(r1-r2))
    return out
def lapla(u):
    return inner(grad(u),grad(u))

def Simulations():

    # nx,ny are number of elements on x and y axis
    # xm, ym are maximum x and y values
    # x0,y0 are minimum x and y values
    # hx, hy are meshwidth on the x and y axis.
    x0=0
    y0=0
    xm=10
    ym=10
    nx=10
    ny=10
    StSt=[0,1]
    perm=0.1

    # FeNiCs variables--------------------------------------------------
    # Initialise the space ---------------------------------------------
    mesh=RectangleMesh(Point(x0,y0),Point(xm,ym),nx,ny)
    #P1 = FiniteElement('P', triangle, 1)
    #element = MixedElement([P1, P1])
    V = FunctionSpace(mesh, 'Lagrange',1)
    # Function(V) indicates a scalar variable over the domain.
    # TestFunction(V) is a scalar test function over the domain for the weak problem.
    u    = Function(V)
    v    = Function(V)
    tmpu0    = Function(V)
    tmpv0    = Function(V)
    q_u      = TestFunction(V)
    q_v      = TestFunction(V)



    #Initialise time step 0
    tmpu0 = Initial(perm,StSt[0],tmpu0)
    tmpv0 = Initial(perm,StSt[1],tmpv0)


    # Define variational problem ---------------------------------------
    F=(lapla(u)+u**2-u*tmpv0)*q_u*dx
    G=(lapla(v)+(tmpu0+v)*(1-v)-v)*q_v*dx

    # Solvers
    problem_u=NonlinearVariationalProblem(F,u, J=derivative(F, u))
    solver_u=NonlinearVariationalSolver(problem_u)
    solver_u.solve()
    tmpu0.assign(u)
    problem_v=NonlinearVariationalProblem(G,v, J=derivative(G, v))
    solver_v=NonlinearVariationalSolver(problem_v)
    solver_v.solve()
    tmpv0.assign(v)
    return

def main(simon=1):
    if simon: Simulations()
    return

if __name__=='__main__':
    main()
    exit()
3

Can you confirm that the line:
problem_u=NonlinearVariationalProblem(F,u, J=derivative(F, u)) is finding u which estimates F=0. I had thought that the second term was to indicate the variable for the solution and the variable being solved for was given by the trial function.
Hence my confusion with two variables.

4

Yep. You can see the source here for example.

The issue with NonlinearVariationalFormulation is that it hides a lot of the actual maths that goes on behind the scenes. The symbolic Gateaux derivative computed by derivative(F, u) implicitly constructs a trial function based on the space in which u is defined yielding a bilinear formulation. You can specify it yourself if you need to as the optional third argument (documentation here).

You may get a better idea of what’s going on by taking a look at this. Or moving to dolfinx which houses the latest developments of dolfin (old dolfin is now very old). dolfinx is far less implicit. In (subjective) order of increasing complexity the following dolfinx demos solve nonlinear problems:

nonlinear Poisson
hyperelasticity
Cahn-Hilliard

5

Thank you, that’s really helpful.

1 Like
6

Just to throw in some more references for DOLFINx to add to @nate’s suggestions: Custom Newton solvers — FEniCSx tutorial

1 Like