Symbolic form of derivated stiffness matrix

1

Hi all,

given is the nonlinear poisson problem

\begin{align*} -\mathrm{div}{(q(x) \nabla u(x))} &= f(x), \quad x \in D = (0,1)^2 \\ u(x) &= 0, \quad x \in \partial D \end{align*}

Let’s use the diffusion coefficient q(x) = 1.0 + r_1 \cdot \sin{(\pi x_1)} + r_2 \cdot \sin{(\pi x_2)} as an example.

The weak formulation is given by

a(u, v) := \int_{D} q \nabla u \nabla v \, \mathrm{d}x = \int_{D} f v \, \mathrm{d}x =: L(v), \quad \forall v \in \hat{V}

I’d like to get the symbolic form of \nabla_r A = (\frac{\partial A}{\partial r_1}, \frac{\partial A}{\partial r_2}) as a python function of r_1 and r_2. Here A denotes the assembled stiffness matrix with A_{ij} = a(\phi_i, \phi_j) for the bilinear form a.

Here’s what I tried so far:

import fenics
import numpy as np

def evaluate_grad_A(r1, r2):
	# rhs of the PDE
	f = fenics.Expression(rhs_string, degree=2)

	# Homogenous Dirichlet Boundary, i.e. u = u_D = 0 on the boundary
	u_D = fenics.Constant(0.0)
	bc = fenics.DirichletBC(V, u_D, lambda x, on_boundary: on_boundary)

	# Define variational problem
	u = fenics.TrialFunction(V)
	v = fenics.TestFunction(V)

	dot, grad, dx = fenics.dot, fenics.grad, fenics.dx
	sin, cos, pi = fenics.sin, fenics.cos, fenics.pi

	x = fenics.SpatialCoordinate(mesh)
	r = fenics.Constant((r1, r2))
	q = fenics.Constant(1.0) + r[0] * sin(pi*x[0]) + r[1] * sin(pi*x[1])

	# Bilinear form
	a = q * dot(grad(u), grad(v)) * dx

	# define r as variable for the derivative
	w = fenics.variable(r)

	dF = fenics.diff(a, w)

	L = f*v * fenics.dx
	u = fenics.Function(V)

	# Assemble
	A, b = fenics.assemble_system(A_form = dF, b_form = L, bcs=bc)
	return A.array()

if __name__ == '__main__':
	fenics.parameters["reorder_dofs_serial"] = False # Disable reordering

	rhs_string = "100*(sin(2*pi*x[1])-sin(2*pi*x[0]))"
	mesh_nx, mesh_ny = 3, 3
	mesh = fenics.UnitSquareMesh(mesh_nx, mesh_ny)
	V = fenics.FunctionSpace(mesh, "P", 1)

Unfortunately, it isn’t working:

In [28]: evaluate_grad_A(0.0, 0.0)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-28-7872e1f48688> in <module>
----> 1 evaluate_grad_A(0.0, 0.0)

<ipython-input-15-88778f7be973> in evaluate_grad_A(r1, r2)
     33
     34         # Assemble
---> 35         A, b = fenics.assemble_system(dF, L, bc)
     36         return A.array()
     37

/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py in assemble_system(A_form, b_form, bcs, x0, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, A_tensor, b_tensor, backend)
    428         assembler.assemble(A_tensor, b_tensor, x0)
    429     else:
--> 430         assembler.assemble(A_tensor, b_tensor)
    431
    432     return A_tensor, b_tensor

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-cZGj9S/dolfin-2019.2.0~git20200629.946dbd3/dolfin/la/PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

Am I missing something? In case you need further information, please let me know. Any help is really appreciated!

Best,
Johnny

2

Note that your code does not run, and requires some modification to be executed. How about the following:

import fenics
import numpy as np
fenics.parameters["reorder_dofs_serial"] = False # Disable reordering

rhs_string = "100*(sin(2*pi*x[1])-sin(2*pi*x[0]))"
mesh_nx, mesh_ny = 3, 3
mesh = fenics.UnitSquareMesh(mesh_nx, mesh_ny)
V = fenics.FunctionSpace(mesh, "P", 1)
# rhs of the PDE
f = fenics.Expression(rhs_string, degree=2)

# Homogenous Dirichlet Boundary, i.e. u = u_D = 0 on the boundary
u_D = fenics.Constant(0.0)
bc = fenics.DirichletBC(V, u_D, lambda x, on_boundary: on_boundary)

# Define variational problem
u = fenics.TrialFunction(V)
v = fenics.TestFunction(V)

dot, grad, dx = fenics.dot, fenics.grad, fenics.dx
sin, cos, pi = fenics.sin, fenics.cos, fenics.pi

x = fenics.SpatialCoordinate(mesh)
r1,r2= fenics.Constant(0),fenics.Constant(0)

q = fenics.Constant(1.0) + r1 * sin(pi*x[0]) + r2 * sin(pi*x[1])

# Bilinear form
a = q * dot(grad(u), grad(v)) * dx

# define r as variable for the derivative

dF1 = fenics.diff(a, r1)
dF2 = fenics.diff(a, r2)

L = f*v * fenics.dx
u = fenics.Function(V)
dF = dF1+dF2
# Assemble
A, b = fenics.assemble_system(A_form = dF, b_form = L, bcs=bc)