Integrating a functional, followed by differentiation

variational formulation
1

Hello,
I am trying to solve the Allen-Cahn equation. My formulation is complex, and I am a bit lost on how to implement it.
First, I have an energy functional E, which has the form:
equation
in which Eta and c_s are both functions, but c_s is a known function, evaluated in the previous time step. For now, it can be considered as a constant function.

Then I have my strong form:
equation(1)
The operator on the right-hand side is a variational differential operator. It can be considered as a differential operator for a functional. So basically, after evaluating E with the first formula, I need to differentiate it with respect to Eta.

My weak form looks like:
equation(2)
with v as the test function, Eta_{n+1} as the trial function and Eta_n is known from previous time step.

My question is: How do I implement this?
First, I need to evaluate the energy functional E over the domain by integration, but maintain Eta as a function. This is where I am lost. I thought of using the assemble() function to evaluate the integral, but this gives me a single scalar. I believe E might be different for each cell and thus can not be a single scalar.
Second, how can I differentiate this functional E with respect to Eta and apply it in the weak form.

I would like to know, how I could construct the variational formulation for my problem and implement it in the solve() function. I have attached my code, but I apologize, I have not been able to make much progress.

Any help is greatly appreciated.

Code:

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

tol = 1e-4

def main():
    T = 150
    num_steps = 10
    dt = T / num_steps

    # Create mesh and define function space
    # mesh = dolfin.Mesh()
    mesh = UnitSquareMesh(8, 8)
   
    V = FunctionSpace(mesh, 'Lagrange', 1)

    eta = Function(V)
    eta = interpolate(Constant(1.0), V)
    c_s = Function(V)
    c_s = interpolate(Constant(1.0), V)

    i = assemble(-2*eta**3*c_s*dx)
    print(i)

if __name__ == '__main__':
    main()
1 Like
2

FEniCS’s UFL includes the derivative(E,eta,v) function, which implements Gateaux differentiation of a functional E with respect to a Function, eta, in the direction v. In mathematical notation, that is

D_vE[\eta] = \left.\frac{d}{d\epsilon}\left(E[\eta + \epsilon v]\right)\right\vert_{\epsilon=0}\text{ ,}

where \epsilon is a scalar. If the direction v is left out (i.e., derivative(E,eta)), then it is taken to be an arbitrary TestFunction in the same FunctionSpace as eta. Based on the Wikipedia article on “functional derivatives”, it seems that the notation from your post is a common alternative way of writing the Gateaux derivative, i.e.,

D_vE = \int_V\frac{\delta E}{\delta \eta} v\,dV\text{ .}
1 Like
3

Thank you very much for your reply. I will try it out and update here.

4

Here is the complete implementation for the above mentioned problem for anybody looking at this later.
Note: This will not give meaningful results or even solve, because I formulated these equations just to ask my question.

from __future__ import print_function
from fenics import *

tol = 1e-14


class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0, tol)


def main():

    # Define simulation parameters
    T = 0.1
    num_steps = 10
    dt = T / num_steps

    # Define mesh
    mesh = UnitSquareMesh(8, 8)

    # Define function space
    V = FunctionSpace(mesh, 'Lagrange', 1)

    # Define dirichlet boundary conditions
    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
    bottom = Bottom()
    bottom.mark(boundaries, 1)
    bc_bottom = DirichletBC(V, Constant(1), boundaries, 1)
    bcs = [bc_bottom]

    # Define initial conditions
    eta_n = project(Constant(1.0), V)
    c_s = project(Constant(1.0), V)

    # Define functions
    eta = Function(V)
    v = TestFunction(V)

    vtkfile1 = File('eta.pvd')

    # Iteration over time
    for i in range(num_steps):

        # Allen-Cahn
        E = -2*eta**3*c_s
        E_derv = derivative(E, eta, v)

        F_eta = -dt*E_derv*dx + eta_n*v*dx
        solve(F_eta == 0, eta, bcs, solver_parameters={"newton_solver": {"relative_tolerance": 1e-6}})

        eta_n.assign(eta)
        vtkfile1 << eta


if __name__ == '__main__':
    main()