Solving coupled PDEs with mixed elements for morpho elasticity

1

Dear all,

I am trying to work out the code for the morpho elasticity below:
\frac{D{\bf v}}{Dt}+{\bf v}(\nabla\cdot{\bf v})-\nabla\cdot{\bf \sigma}={\bf f},
{\bf \sigma}=E{\bf \epsilon},
\frac{D{\bf \epsilon}}{Dt}+{\bf \epsilon}skw(\nabla{\bf v})-skw(\nabla{\bf v}){\bf \epsilon}+(tr({\bf \epsilon})-1)sym(\nabla{\bf v})={\bf 0}.

Here, we are looking for the solution to ({\bf \epsilon,v}) with homogeneous boundary condition for {\bf v}, where {\bf \epsilon} is a tensor (or 2*2 matrix in 2D) and {\bf v} is a vector for velocity.

By using Reynold’s transport theorem and for P1 elements, I got the weak form and write the code like

from future import print_function

import necessary packages

import numpy as np
from array import *
from pylab import *

import FEniCS packages

from fenics import *
from mshr import *
from dolfin import *

##########################################################

Parameters’ Value

##########################################################
Es=100
##########################################################

settings of cell moving

##########################################################

computational domain

x0=60
y0=40
CompDomain=Rectangle(Point(-x0, -y0), Point(x0, y0))
mesh = generate_mesh(CompDomain, 40)

meshpoints=mesh.coordinates()

parameters[‘reorder_dofs_serial’] = False

v_1=VectorElement(‘BDM’,triangle, 1)
v_2=VectorElement(‘P’,triangle, 1)
v=MixedElement([v_1,v_2])
V=FunctionSpace(mesh, v)

v_D = Constant((0,0))
f=Constant((1,0))

def boundary(x, on_boundary):
return on_boundary

bc = DirichletBC(V.sub(1), v_D, boundary)

def sigma(v):
return Es*eps

Define trial and test functions

(eps,v) = TrialFunctions(V)
(tau, phi) = TestFunctions(V)

F=inner(sigma(v),grad(phi))*dx+inner(v,phi)*ds-inner(f,phi)dx+inner(epsskew(grad(v)),tau)*dx-inner(skew(grad(v))*eps,tau)*dx+(tr(eps)-1)*inner(sym(grad(v)),tau)*dx
w=Function(V)
solve(F==0, w,bc)

Then with current code I got the error:
Traceback (most recent call last):
File “/usr/lib/python3/dist-packages/dolfin/compilemodules/jit.py”, line 142, in jit
result = ffc.jit(ufl_object, parameters=p)
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 218, in jit
module = jit_build(ufl_object, module_name, parameters)
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 134, in jit_build
generate=jit_generate)
File “/usr/lib/python3/dist-packages/dijitso/jit.py”, line 167, in jit
header, source, dependencies = generate(jitable, name, signature, params[“generator”])
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 67, in jit_generate
prefix=module_name, parameters=parameters, jit=True)
File “/usr/lib/python3/dist-packages/ffc/compiler.py”, line 143, in compile_form
prefix, parameters, jit)
File “/usr/lib/python3/dist-packages/ffc/compiler.py”, line 185, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
File “/usr/lib/python3/dist-packages/ffc/analysis.py”, line 94, in analyze_ufl_objects
for form in forms)
File “/usr/lib/python3/dist-packages/ffc/analysis.py”, line 94, in
for form in forms)
File “/usr/lib/python3/dist-packages/ffc/analysis.py”, line 178, in _analyze_form
do_apply_restrictions=True)
File “/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py”, line 388, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is
File “/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py”, line 152, in check_form_arity
check_integrand_arity(itg.integrand(), arguments)
File “/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py”, line 145, in check_integrand_arity
args = map_expr_dag(rules, expr, compress=False)
File “/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py”, line 37, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File “/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py”, line 86, in map_expr_dags
r = handlers[v.ufl_typecode](v, *[vcache[u] for u in v.ufl_operands])
File “/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py”, line 42, in sum
raise ArityMismatch(“Adding expressions with non-matching form arguments {0} vs {1}.”.format(a, b))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement(‘Lagrange’, triangle, 1), dim=2), 3), MixedElement(VectorElement(FiniteElement(‘Brezzi-Douglas-Marini’, triangle, 1), dim=2), VectorElement(FiniteElement(‘Lagrange’, triangle, 1), dim=2))), 1, None),).
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File “/home/qiyaopeng/Downloads/pycharm-community-2017.3.3/helpers/pydev/pydev_run_in_console.py”, line 53, in run_file
pydev_imports.execfile(file, globals, locals) # execute the script
File “/home/qiyaopeng/Downloads/pycharm-community-2017.3.3/helpers/pydev/_pydev_imps/_pydev_execfile.py”, line 18, in execfile
exec(compile(contents+“\n”, file, ‘exec’), glob, loc)
File “/home/qiyaopeng/Desktop/PhD/StudyNotes/Project/ProjectCode/Morpho_Elas_Approach/morpho_2D.py”, line 90, in
solve(F==0, w,bc)
File “/usr/lib/python3/dist-packages/dolfin/fem/solving.py”, line 300, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/lib/python3/dist-packages/dolfin/fem/solving.py”, line 344, in _solve_varproblem
form_compiler_parameters=form_compiler_parameters)
File “/usr/lib/python3/dist-packages/dolfin/fem/solving.py”, line 149, in init
F = Form(F, form_compiler_parameters=form_compiler_parameters)
File “/usr/lib/python3/dist-packages/dolfin/fem/form.py”, line 92, in init
mpi_comm=mesh.mpi_comm())
File “/usr/lib/python3/dist-packages/dolfin/compilemodules/jit.py”, line 70, in mpi_jit
return local_jit(*args, **kwargs)
File “/usr/lib/python3/dist-packages/dolfin/compilemodules/jit.py”, line 147, in jit
“ffc.jit failed with message:\n%s” % (tb_text,))
File “/usr/lib/python3/dist-packages/dolfin/cpp/common.py”, line 2808, in dolfin_error
return _common.dolfin_error(location, task, reason)
RuntimeError:
*** -------------------------------------------------------------------------
*** 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 perform just-in-time compilation of form.
*** Reason: ffc.jit failed with message:
Traceback (most recent call last):
File “/usr/lib/python3/dist-packages/dolfin/compilemodules/jit.py”, line 142, in jit
result = ffc.jit(ufl_object, parameters=p)
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 218, in jit
module = jit_build(ufl_object, module_name, parameters)
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 134, in jit_build
generate=jit_generate)
File “/usr/lib/python3/dist-packages/dijitso/jit.py”, line 167, in jit
header, source, dependencies = generate(jitable, name, signature, params[“generator”])
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 67, in jit_generate
prefix=module_name, parameters=parameters, jit=True)
File “/usr/lib/python3/dist-packages/ffc/compiler.py”, line 143, in compile_form
prefix, parameters, jit)
File “/usr/lib/python3/dist-packages/ffc/compiler.py”, line 185, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
File “/usr/lib/python3/dist-packages/ffc/analysis.py”, line 94, in analyze_ufl_objects
for form in forms)
File “/usr/lib/python3/dist-packages/ffc/analysis.py”, line 94, in
for form in forms)
File “/usr/lib/python3/dist-packages/ffc/analysis.py”, line 178, in _analyze_form
do_apply_restrictions=True)
File “/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py”, line 388, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is
File “/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py”, line 152, in check_form_arity
check_integrand_arity(itg.integrand(), arguments)
File “/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py”, line 145, in check_integrand_arity
args = map_expr_dag(rules, expr, compress=False)
File “/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py”, line 37, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File “/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py”, line 86, in map_expr_dags
r = handlers[v.ufl_typecode](v, *[vcache[u] for u in v.ufl_operands])
File “/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py”, line 42, in sum
raise ArityMismatch(“Adding expressions with non-matching form arguments {0} vs {1}.”.format(a, b))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement(‘Lagrange’, triangle, 1), dim=2), 3), MixedElement(VectorElement(FiniteElement(‘Brezzi-Douglas-Marini’, triangle, 1), dim=2), VectorElement(FiniteElement(‘Lagrange’, triangle, 1), dim=2))), 1, None),).
.
*** Where: This error was encountered inside jit.py.
*** Process: 0

*** DOLFIN version: 2017.2.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

If I take v_1 as P1 element, the shape of inner product does not match. Can anyone help me out of this? Many thanks!

Best regards,
Alice

2

The problem here is that you are using a TrialFunction for the solution to a nonlinear problem. TrialFunction should only be used in linear problems. You should use a Function for the solution when defining the form for a nonlinear residual. (But, once I fix that, the nonlinear residual is nan on the second Newton iteration, so there are still other problems to debug.)

3

Hi Kamensky,

Thanks for your comment and I also get the error but in the 0 iteration already:

Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 2.033e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
File “/home/qiyaopeng/Downloads/pycharm-community-2017.3.3/helpers/pydev/pydev_run_in_console.py”, line 53, in run_file
pydev_imports.execfile(file, globals, locals) # execute the script
File “/home/qiyaopeng/Downloads/pycharm-community-2017.3.3/helpers/pydev/_pydev_imps/_pydev_execfile.py”, line 18, in execfile
exec(compile(contents+"\n", file, ‘exec’), glob, loc)
File “/home/qiyaopeng/Desktop/PhD/StudyNotes/Project/ProjectCode/Morpho_Elas_Approach/morpho_2D.py”, line 91, in
solve(F==0, w,bc)
File “/usr/lib/python3/dist-packages/dolfin/fem/solving.py”, line 300, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/lib/python3/dist-packages/dolfin/fem/solving.py”, line 349, in _solve_varproblem
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** 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 successfully call PETSc function ‘MatSetValuesLocal’.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /build/dolfin-fyfCiA/dolfin-2017.2.0.post0/dolfin/la/PETScMatrix.cpp.
*** Process: 0

*** DOLFIN version: 2017.2.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

Does that mean I should try other iteration method or something else? Thanks in advance.

KInd regards,
Alice

4

The reason you’re getting a different error is probably because you didn’t remove the line

w = Function(V)

right before the solve. You want the Function passed as the solution argument to solve to be the same Function used as the solution when defining the residual.

5

Ah I see. Thanks a lot for that and sorry that I never worked with nonlinear problem before, I should have checked the tutorial more carefully.

Besides, I also noticed that I got nan since the second iteration, which might caused by wrong weak form.