Error in Changing Constant Coefficient to 3x3 Matrix Coefficient

1

I am attempting to solve a nonlinear version of Maxwell’s equations and was able to make the following minimal working example code work. However, instead of keeping the material coefficient ‘eps’ constant the way it is in code, I now want to make it a 3x3 coefficient matrix. My attempt at defining the 3x3 coefficient matrix is shown in the comments. However, this results in FEniCS showing an error when it gets to my weak formulation ‘F’. These errors are exemplified after the code below:

from dolfin import *
import numpy as np

T = 5.0            # final time
num_steps = 500    # number of time steps
dt = T / num_steps # time step size
k = Constant(dt)

# Material Coefficients (SCALAR)
eps = Constant(0.1) # GOAL TO MAKE THIS A 3x3 Tensor using: Constant(np.array([[0.98, 0, 0], [0, 0.19, 0], [0, 0, 0.19]])) 

mesh = UnitCubeMesh(1,1,1)
P1 = FiniteElement("N1curl", mesh.ufl_cell(), 1) 
TH = P1*P1
V = FunctionSpace(mesh, TH)

# Define test functions
test_E, test_H = TestFunctions(V)

# Define functions for velocity and concentrations
u = Function(V)
u_n = Function(V)

# Split system functions to access components
trial_E, trial_H = split(u)
prev_E, prev_H = split(u_n)

F = eps*inner(trial_E - prev_E, test_E)/k*dx 

# Create VTK files for visualizing output
vtkfile_u_1 = File('reaction_system/u_5.pvd')
vtkfile_u_2 = File('reaction_system/u_6.pvd')


# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Solve variational problem for time step
    solve(F == 0, u) 

    # Save solution to file (VTK)
    _u_1, _u_2 = u.split()
    vtkfile_u_1 << (_u_1, t)
    vtkfile_u_2 << (_u_2, t)

    # Update previous solution
    u_n.assign(u)
    plot(_u_1)

Errors
Following are the various errors I get for replacing the definition of F in the code above. It is the same code all throughout, except for the line with F in it.

  1. When I define F the way it is in the code above (but with eps defined as a 3x3 matrix), I get the error:

    Can only integrate scalar expressions. 
The integrand is a tensor expression with value shape (3, 3) and free indices with labels ().

  2. If I replace F with the following (inner product):

    F = inner(eps,inner(trial_E - prev_E, test_E))/k*dx

    I get

    Shapes do not match: <Coefficient id=4934954504> and <Inner id=4935313448>.

  3. I get a similar error as 2 when I switch the order of the product to:

    F = inner(inner(trial_E - prev_E, test_E),eps)/k*dx

  4. If I try a dot product:

    F = dot(inner(trial_E - prev_E, test_E),eps)/k*dx

    I get:

    Dot product requires non-scalar arguments, got arguments with ranks 0 and 2.

So, either the way I defined the 3x3 matrix for ‘eps’ (shown in comments in code above) is incorrect, or the way I defined F.
Any advice? Thanks!

2

I found a very similar error here: /qa/7346/an-error-in-assembling-a-variational-problem/

So, things I learned are:

  1. for inner(), the objects should be the same size. I should just use:

    F=eps*inner(trial_E - prev_E, test_E)/k*dx

  2. But this gives me the expected error. So, the problem is that I should specify my goal better.

My hope was to specify a 3x3 material coefficient matrix. According to the given answer, I should then handle the system component-wise.

  1. But this means that a simple fix would be to move eps inside the inner product!

    F=inner(eps*(trial_E - prev_E), test_E)/k*dx

    However, I am not sure about how this changes my system. Any clarification would be appreciated!

Based off this, for simplicity, I will just use a vector coefficient and move on. However, if anyone has ideas to handle this problem effectively (just for curiosity), do let me know! Thanks!

3

You can do the following for eps:

eps = Constant(((0.98, 0.0, 0.0), (0.0, 0.19, 0.0), (0.0, 0.0, 0.19)))

and use this in the weak form as

F = inner(eps*(trial_E - prev_E), test_E)/k*dx

which is what you found out already. Here, it does not make a difference, since the matrix is symmetric, for non-symmetric ones you have to be careful which side you multiply by it.

For example, take a look at the fenics book (chapter 17.4) for the implementation of such routines