Derive weak form from variational principle

variational formulation
1

Hi all,

I’m trying to solve a set of coupled non-linear PDEs in FEniCSx, but I’m having trouble understanding how to write the variational problem. For now I’m assuming Dirichlet BCs, although later I will need to include rather complicated Robin BCs. Moreover, the problem is such that spherical coordinates and basis are the natural set to work on.

First of all, we have an energy functional F that depends on the vector order parameter \mathbf{m}\equiv(m_r,m_{\theta},m_{\phi}) and on its derivatives, which for simplicity I note as \partial\mathbf{m}:

F[\mathbf{m},\partial\mathbf{m}]= \int_{\Omega} \mathrm{d}V f_B(\mathbf{m},\partial\mathbf{m}) = \int_{\Omega} \mathrm{d}r \mathrm{d}\theta \mathrm{d}\phi\, J(r, \theta) f_B(\mathbf{m},\partial\mathbf{m})

with J(r, \theta)=r^2\sin\theta. Without specifying f_B, I can impose \delta F =0 to minimise the energy, and thus derive the Euler-Lagrange equations for m_i,\ i=r,\theta,\phi:

J(r, \theta)\left[\dfrac{\partial f_B}{\partial m_i}-\nabla\cdot \dfrac{\partial f_B}{\partial (\partial m_i)}\right]=0

where \nabla\cdot is the divergence in spherical coordinates and

\dfrac{\partial f_B}{\partial (\partial m_i)}\equiv \left( \dfrac{\partial f_B}{\partial (\partial_r m_i)}, \dfrac{\partial f_B}{\partial (\partial_{\theta} m_i)},\dfrac {\partial f_B}{\partial (\partial_{\phi} m_i)}\right)

is a vector containing the derivatives of f_B with respect to the derivatives of m_i.

******First question: can I derive the weak formulation directly from this general equation or (even better) from the energy functional, and then replace with f_B? Or is it necessary to specify f_B from the beginning?

In this case,

f_B=\dfrac{a}{2}m^2+\dfrac{c}{4}m^4 + d \left\{\left(\nabla\cdot\mathbf m\right)^2 + \left[\mathbf m\cdot \left(\nabla\times\mathbf m\right)\right]^2 + \left[\mathbf m\times \left(\nabla\times\mathbf m\right)\right]^2\right\} \,.

Despite assuming rotational symmetry (so nothing depends on \phi) and taking m_{\phi}=0, this gives rise to very complicated equations for the remaining components of \mathbf{m}\equiv(m_r,m_{\theta},0):

image
image940×289 10.7 KB

(As an image from Mathematica so I don’t make any mistakes while copying them.)

******Second question: Any suggestions on how to derive the weak form efficiently and then implement it in FEniCSx?

Thanks in advance!

2

The canonical approach in FEniCS would be to specify F in UFL, then use the derivative() function to automatically take its Gateaux derivative with respect to \mathbf{m} in the direction of a test function \delta\mathbf{m} to get \delta F, which would then be the residual of the weak problem. This is done in the Legacy FEniCS hyperelasticity demo, which should be easy to adapt to FEniCSx. (The current FEniCSx hyperelasticity demo formulates the problem in a different but equivalent way, by taking a partial derivative of the energy density with diff(), to get the first Piola–Kirchhoff stress; this is slightly more work, but can be clearer when combining with time integration and non-conservative forcing. That would be analogous to starting from your Euler–Lagrange form, integrating by parts, and only automating the partials of f_B with respect to \mathbf{m} and \nabla\mathbf{m}.)

To handle the change of coordinates, it would probably be cleanest to define auxiliary Python functions for the polar divergence and curl, interpreting UFL’s SpatialCoordinate to be (r,\theta,\phi), where UFL grad() would be a vector of partial derivatives with respect to polar coordinates (and not a true coordinate-independent gradient) and a dx-type Measure would similarly be dr\,d\theta\,d\phi (instead of a physical volume measure).

1 Like
3

Hi kamensky, thanks for your reply. I’ve started by defining functions for the spherical curl, divergence and gradient, as you suggested, but I don’t know how to return a UFL type instead of a list or numpy array. As a consequence, I can’t compute dot products between a vector function and its curl, for example. Here’s a mwe:

import ufl
from dolfinx import mesh, fem
from mpi4py import MPI

dim = 3 # dimension of the OP

domain = mesh.create_unit_cube(MPI.COMM_WORLD,10,10,10)
x = ufl.SpatialCoordinate(domain)
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, dim)
V = fem.FunctionSpace(domain, element) # space of vector functions with dim components

m = ufl.TrialFunction(V) # dim-dimensional vector OP

def curl_sph(A, x):
    xr, xq, xf = x; jr, jq, jf = 0, 1, 2
    Ar, Aq, Af = A
    curl_r = (ufl.Dx(Af*ufl.sin(xq),jq)-ufl.Dx(Aq,jf))/(xr*ufl.sin(xq))
    curl_q = (ufl.Dx(Ar,jf)/ufl.sin(xq)-ufl.Dx(xr*Af,jr))/xr
    curl_f = (ufl.Dx(xr*Aq,jr)-ufl.Dx(Ar,jq))/xr
    return [curl_r, curl_q, curl_f]

#ufl.dot([1,0,0],[0,1,0])
ufl.dot(m,curl_sph(m,x))

---------------------------------------------------------------------------
UFLValueError                             Traceback (most recent call last)
Input In [1], in <cell line: 23>()
     20     return [curl_r, curl_q, curl_f]
     22 #ufl.dot([1,0,0],[0,1,0])
---> 23 ufl.dot(m,curl_sph(m,x))

File /usr/local/lib/python3.10/dist-packages/ufl/operators.py:173, in dot(a, b)
    171 "UFL operator: Take the dot product of *a* and *b*. This won't take the complex conjugate of the second argument."
    172 a = as_ufl(a)
--> 173 b = as_ufl(b)
    174 if a.ufl_shape == () and b.ufl_shape == ():
    175     return a * b

File /usr/local/lib/python3.10/dist-packages/ufl/constantvalue.py:437, in as_ufl(expression)
    435     return IntValue(expression)
    436 else:
--> 437     raise UFLValueError("Invalid type conversion: %s can not be converted"
    438                         " to any UFL type." % str(expression))

UFLValueError: Invalid type conversion: [Division(Sum(Indexed(Grad(Product(Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(2),))), Sin(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(1),)))))), MultiIndex((FixedIndex(1),))), Product(IntValue(-1), Indexed(Grad(Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(1),)))), MultiIndex((FixedIndex(2),))))), Product(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(0),))), Sin(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(1),)))))), Division(Sum(Product(IntValue(-1), Indexed(Grad(Product(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(0),))), Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(2),))))), MultiIndex((FixedIndex(0),)))), Division(Indexed(Grad(Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(0),)))), MultiIndex((FixedIndex(2),))), Sin(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(1),)))))), Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(0),)))), Division(Sum(Indexed(Grad(Product(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(0),))), Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(1),))))), MultiIndex((FixedIndex(0),))), Product(IntValue(-1), Indexed(Grad(Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(0),)))), MultiIndex((FixedIndex(1),))))), Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(0),))))] can not be converted to any UFL type.

In fact, this extremely minimal example raises the same error:

import ufl
ufl.dot([1,0,0],[0,1,0])

---------------------------------------------------------------------------
UFLValueError                             Traceback (most recent call last)
Input In [1], in <cell line: 3>()
      1 import ufl
----> 3 ufl.dot([1,0,0],[0,1,0])

File /usr/local/lib/python3.10/dist-packages/ufl/operators.py:172, in dot(a, b)
    170 def dot(a, b):
    171     "UFL operator: Take the dot product of *a* and *b*. This won't take the complex conjugate of the second argument."
--> 172     a = as_ufl(a)
    173     b = as_ufl(b)
    174     if a.ufl_shape == () and b.ufl_shape == ():

File /usr/local/lib/python3.10/dist-packages/ufl/constantvalue.py:437, in as_ufl(expression)
    435     return IntValue(expression)
    436 else:
--> 437     raise UFLValueError("Invalid type conversion: %s can not be converted"
    438                         " to any UFL type." % str(expression))

UFLValueError: Invalid type conversion: [1, 0, 0] can not be converted to any UFL type.

How can I return something compatible with ufl.dot?

4

Use ufl.as_vector([curl_r, curl_q, curl_f]).

1 Like
5

Hi again. So thanks to @dokken I’ve been able to write the energy, but now I can’t take the Gateaux derivative as suggested by @kamensky. I’ve followed the steps in Hyperelasticity. Here’s a MWE:

import ufl
from dolfinx import mesh, fem
from mpi4py import MPI

dim = 3 # dimension of the OP
domain = mesh.create_unit_cube(MPI.COMM_WORLD,10,10,10)
x = ufl.SpatialCoordinate(domain)
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, dim)
V = fem.FunctionSpace(domain, element) # space of vector functions with dim components

def Jacob(x):
    """Jacobian of the transformation from Cartesian to spherical coordinates."""
    xr, xq, xf = x
    return xr**2*ufl.sin(xq)

m = ufl.TrialFunction(V) # dim-dimensional vector OP
phi = ufl.TestFunction(V) # test function

F = ufl.dot(m,m)*Jacob(x)*ufl.dx 
dF = ufl.derivative(F, m, phi) # variation of the energy

which returns the error:

ERROR:UFL:Invalid coefficient type for <Argument id=140129317279424>
---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
Input In [2], in <cell line: 23>()
     20 phi = ufl.TestFunction(V) # test function
     22 F = ufl.dot(m,m)*Jacob(x)*ufl.dx 
---> 23 dF = ufl.derivative(F, m, phi)

File /usr/local/lib/python3.10/dist-packages/ufl/formoperators.py:272, in derivative(form, coefficient, argument, coefficient_derivatives)
    248 def derivative(form, coefficient, argument=None, coefficient_derivatives=None):
    249     """UFL form operator:
    250     Compute the Gateaux derivative of *form* w.r.t. *coefficient* in direction
    251     of *argument*.
   (...)
    269     ``Coefficient`` instances to their derivatives w.r.t. *coefficient*.
    270     """
--> 272     coefficients, arguments = _handle_derivative_arguments(form, coefficient,
    273                                                            argument)
    275     if coefficient_derivatives is None:
    276         coefficient_derivatives = ExprMapping()

File /usr/local/lib/python3.10/dist-packages/ufl/formoperators.py:222, in _handle_derivative_arguments(form, coefficient, argument)
    220 else:
    221     if not isinstance(c, Indexed):
--> 222         error("Invalid coefficient type for %s" % ufl_err_str(c))
    223     f, i = c.ufl_operands
    224     if not isinstance(f, Coefficient):

File /usr/local/lib/python3.10/dist-packages/ufl/log.py:158, in Logger.error(self, *message)
    156 "Write error message and raise an exception."
    157 self._log.error(*message)
--> 158 raise self._exception_type(self._format_raw(*message))

UFLException: Invalid coefficient type for <Argument id=140129317279424>

Does anyone know what the problem with phi thus defined is?

6

m should probably be a dolfinx.fem.Function

7

Thanks @dokken for your help! I have more problems now, but I think it’s better to start a new topic.