How to write Coupled, non-linear gradients in variational formula

variational formulation
1

Hi, everyone. I’m trying to write my variational formula that has coupled terms that are derivatives of scalars, u and phi2,. The terms are:
u(x,t) \nabla (\nabla \phi_2(x,t))
\nabla u(x,t) \nabla phi_2(x,t)
I’ve tried formulating them several different ways, and I had issues with taking non-scalars in my dot command, so I dot the gradient with a FacetNormal function, N, but I get:
“ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal.”
Here’s my code snipped containing these terms in variational formula

# Create mesh and functions
mesh = IntervalMesh(100,0,1)
V = FunctionSpace(mesh,'P',1)
Q = FunctionSpace(mesh,'P',1)
phi2 = Function(Q)
w = TestFunction(Q)
u = Function(V)
v = TestFunction(V)

# Create Initial values
# Define initial values
phi2_0 = Constant(phi0)
phi2_n = project(phi2_0, Q) 
u_0 = Constant(0)
u_n = project(u_0, V) 

# Define normal unit
N = FacetNormal(mesh)

F1 = 	dot(N*u,grad( dot(grad(phi2_n), N) ))*v*dx \
	dot(grad(u),grad(phi2_n))*v*dx
F2 = 	 dot( N*u_n, grad(dot(grad(phi2), N)) )*w*dx \
	+ dot( grad(u_n),grad( dot(grad(phi2), N) ) )*w*dx

I’ve removed the lines from F1 and F2 that are not causing problems and not included the extra code for the BC’s and constants. My info for FEniCS: Ubuntu 20.04 LTS (GNU/Linux 4.4.0-19041-Microsoft x86_64). I realize usually I can integrate by parts but I don’t see how that works for two coupled, non-linear functions, especially when they have derivatives.

2

The reason for having two F equations is that there are two PDE’s for these functions, but this coupled term shows up in both.

3

You cannot use facetnormals for cell integrals, as they are not defined. Note that you should integrate your equation by parts, to avoid having grad^2 terms

4

Ah, I did not know not to use facetnormals in integrals. And for getting rid of grad^2 terms, if I integrate that term by parts, then wouldn’t I get
- \nabla u \nabla \phi_2 \nabla v
? In which case would cancel out the other term I need to include?

5

No, as you are using different test functions that wouldn’t cancel out. I would like to note that if you are solving a coupled problem, you should use a mixed function space. Please carefully write your problem by hand with different test functions for each equation and integrate by parts

6

If both of my PDE’s contain the terms
u \nabla (\nabla \phi_2) and \nabla u \nabla \phi_2
wouldn’t they still cancel? Multiplying both these terms by the u test function will cancel them in the equation for u, and then multiplying the two terms in the other equation by the phi2 test function would cancel them, even if they’re added together.
Situations where I see that not happening are

  1. The integration by parts in the u test function can somehow produce the phi2 test function
  2. I combine both PDE’s into one equation and then multiply every term by both test functions.
    Please let me know if I’m missing something obvious.
7

As you have not supplied the strong form of your PDE, I cannot say anything about how you should integrate by parts

8

Green’s first and second identities are probably all you need for standard elliptic problems.

1 Like
9

That looks very helpful, thank you! For reference, here is my set of PDE’s:
\partial u / \partial t = p* \nabla (\nabla u) + q*\nabla (u * \nabla \phi_2)
\nabla (u * \nabla \phi_2) = -r*u*(s - \phi_2)
where p, q, r, and s are constants.

10

For the Green’s identities, it looks like it’s integrating the terms with respect to dV without the multiplication of another function. Is it still valid since my equations will be integrating the two terms both multiplied by the test function, v?
My situation is \int_\Omega \nabla \cdot (u \nabla \phi_2 )v dx
instead of \int_\Omega \nabla \cdot (u \nabla \phi_2 )dx

11

See for instance greens second identity, where there are three variables present. These identities should help you. I would also refer you to having a look at the Stokes equations such that you understand how to deal with mixed problems.

1 Like
12

That example is very helpful, thanks! It seems like the only difference would be that I should define u and \phi_2 as Functions instead of TrialFunctions since they’re nonlinear.
For writing the variational problem, I’ve messed around with a few rearrangements of those two identities. It seems like if I want to get rid of the second derivative for Equation 1
\int_\Omega \nabla \cdot (u \nabla \phi_2 )v dx = \int_\Omega (u \nabla \cdot (v \nabla \phi_2 ) + v \nabla u \cdot \nabla \phi_2 ) dx
then using the logic from Theorem 2
\int_\Omega (u \nabla \cdot (v \nabla \phi_2 ) - \phi_2 \nabla \cdot (v \nabla u ) )dx = \int_{d\Omega} v(u \partial \phi_2 / \partial n - \phi_2 \partial u / \partial n ) ds
it would seem the degree of derivative decreases when multipled by the test function and integrated. Thus the first term of Equation 1 becomes
\int_{d\Omega} v(u * \partial \phi_2 / \partial n) ds
and the second term becomes
\int_{d\Omega} v(u *\partial \phi_2 / \partial n+ u * \partial \phi_2 / \partial n ) ds
Now I have gotten rid of second derivatives. Then in summary:
\int_\Omega \nabla \cdot (u \nabla \phi_2 )v dx = \int_{d\Omega} v(2* u *\partial \phi_2 / \partial n+ u * \partial \phi_2 / \partial n ) ds
This route sounds strange to me, but I combined the two theorems several ways and couldn’t figure the reasonable thing to do with the test function multiplier without maintaining the double derivative.

13

See for instance page 3 of https://folk.uio.no/kent-and/sommerskole/material/book_20april.pdf Where kappa takes the place if u in your calculation.
As i have said multiple times now, you need to use a unique test function for each problem, similar to what is done in the stokes problem.

14

I see; I wasn’t sure if that integration holds if kappa is a function instead. Thank you.
I should’ve written out my entire equation then to clarify how I believe to write it. For the 1-spatial dimension PDE system with p, q, r, s as constants:
\partial u / \partial t = p \nabla \cdot (\nabla u) + q \nabla (u \nabla \phi_2 )
\nabla (u \nabla \phi_2 ) = -r u (s - \phi_2 )
I multiply the first equation by test function v and second equation by test function w, then integrate, resulting in
\int_\Omega ( (u-u_n)v /dt + p \nabla u \cdot \nabla v+ q u \nabla \phi_2 \cdot \nabla v) dx - q \int_{d\Omega} u (\partial \phi_2 /\partial n)v ds \\ - \int_\Omega ( u \nabla \phi_2 \cdot \nabla w + r u (s - \phi_2 )w dx) + \int_{d\Omega} u (\partial \phi_2 /\partial n)w ds
And then the subdomain {d\Omega} terms go to 0 as it usually does because the test functions are 0 at the boundaries. If not, I’m also not sure the difference between \nabla \phi_2 and \partial \phi_2 / \partial n for a 1D problem. The functions would be defined as such:

# Make Mixed space with previously defined distance, L
mesh = IntervalMesh(100,0,L)
P1 = FiniteElement('P', mesh.ufl_cell(), 1) 	
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
# Define test and regular functions
o = Function(V)
u, phi2 = split(o)
v, w = TestFunctions(V)

with p, q, r, and s as constants. I’m working on the code now to write u_n so I can have initial conditions for \phi_2 and u.
Edit made to explain L in mesh and correct integral.

15

I found an example that shares one PDE with me:
https://fenicsproject.org/qa/13269/scalar-mixed-function-space-problem/
I’ve updated its code to use current FEniCS terms, and it ran fine in 2D. However, when I switch into 1D, the Newton Solver won’t converge. Error message:

Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset:  unknown

My code is:

from dolfin import *
import mshr


# Define boundaries; made 1D, getting rid of y-direction
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)


class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)


# Initialize sub-domain instances
left = Left()
right = Right()

# Mesh. Not sure on this number of cells; have tried up to 10,000 and low as 20.
mesh = IntervalMesh(100, 0, 1)

# Marking boundary domains. 
# Got rid of y-boundaries and re-numbered left & right boundaries.
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)

CG1_elem = FiniteElement("CG", mesh.ufl_cell(), 1)

#Defining the function spaces for the concentration
V_c = FunctionSpace(mesh, CG1_elem)

#Defining the function spaces for the voltage
V_phi = FunctionSpace(mesh, CG1_elem)

#Defining the mixed function space
W_elem = MixedElement([CG1_elem, CG1_elem])
W = FunctionSpace(mesh, W_elem)

#Defining the "Trial" functions
u = Function(W)
c, phi = split(u)

#Defining the test  functions
(v_1, v_2) = TestFunctions(W)

# Time variables
dt = 0.05; t = 0; T = 10

# Previous solution
W_const = Constant(1.0)
C_previous= interpolate(W_const, V_c)

# Define Dirichlet boundary conditions at right boundary,
# previously top boundary; (V=4.0 volt)
Voltage = Constant(4.0)
bc_right = DirichletBC(W.sub(0), Voltage, boundaries, 2)

# Define Dirichlet boundary conditions at left boundary,
# previously bottom boundary (V=0.0 volt)
bc_left = DirichletBC(W.sub(0), 0.0, boundaries, 1)

#Collecting the boundary conditions
bcs = [bc_right, bc_left]

#Define eps in the weak form
eps=0.0082

# Define variational form
a = dot(grad(phi), grad(v_2)) * dx \
    + c*v_1 * dx \
    - (1./(2 * (pow(eps,2)))) * c * v_2 * dx \
    + dt * eps * dot(grad(c),grad(v_1)) * dx \
    + dt * eps * c * dot(grad(phi),grad(v_1)) * dx

L = C_previous * v_1 * dx - (1./(2 * (pow(eps,2)))) * v_2 * dx

# For nonlinear solver
F = a - L

# Create VTK files for visualization output
vtkfile_phi = File('solution/phi.pvd')
vtkfile_c = File('solution/c.pvd')


while t <= T:
    solve(F == 0, u, bcs)

    (c, phi) = u.split(True)
    C_previous.assign(c)
    
    # Save solution to file (VTK)
    vtkfile_c << (c, t)
    vtkfile_phi << (phi, t)


    # Update current time
    t += dt

Is there a reason the one spatial dimension version can’t work? I’ve tried using meshes with fewer points and many more points (20 or 1,000 or even 10,000).