ODE Formulation doesn't work for mixed-implicit-explicit time integration schemes

[Finbar_Argus]

I am trying to solve an ODE of the form

\dot{x} = Ax + Bu

I discritise in time to get the following
(x - x_n)/dt = A\hat{x} + Bu

Where \hat{x} depends on the integration scheme. For Explicit Euler \hat{x} = x_n and for Implicit Euler \hat{x} = x.

I want to use a scheme that has \hat{x} = [x_0, x_1, x_{n2}]. I.e a combination of explicit variables and implicit variables. (Symplectic Euler Scheme). However when I try to implement this I receive the following error…

Traceback (most recent call last):
  File "/home/finbar/Documents/git_projects/portHamiltonian_FEM/fenics_projects/wave_eqn/electroMechanical_forDiscourse.py", line 114, in <module>
    a = lhs(F)
  File "/usr/lib/python3/dist-packages/ufl/formoperators.py", line 79, in lhs
    return compute_form_lhs(form)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/formtransformations.py", line 379, in compute_form_lhs
    return compute_form_with_arity(form, 2)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/formtransformations.py", line 345, in compute_form_with_arity
    return map_integrands(_transform, form)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 39, in map_integrands
    for itg in form.integrals()]
  File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 39, in <listcomp>
    for itg in form.integrals()]
  File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 46, in map_integrands
    return itg.reconstruct(function(itg.integrand()))
  File "/usr/lib/python3/dist-packages/ufl/algorithms/formtransformations.py", line 341, in _transform
    e, provides = pe.visit(e)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 114, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/formtransformations.py", line 145, in sum
    part, term_provides = self.visit(term)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 114, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/formtransformations.py", line 145, in sum
    part, term_provides = self.visit(term)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 114, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/formtransformations.py", line 275, in linear_indexed_type
    part, provides = self.visit(expression)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 110, in visit
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 110, in <listcomp>
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 114, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/formtransformations.py", line 275, in linear_indexed_type
    part, provides = self.visit(expression)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 114, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/formtransformations.py", line 275, in linear_indexed_type
    part, provides = self.visit(expression)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 114, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/formtransformations.py", line 275, in linear_indexed_type
    part, provides = self.visit(expression)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 110, in visit
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 110, in <listcomp>
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 114, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/formtransformations.py", line 275, in linear_indexed_type
    part, provides = self.visit(expression)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py", line 110, in visit
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/lib/python3/dist-packages/ufl/algorithms/formtransformations.py", line 309, in list_tensor
    error("PartExtracter does not know how to handle list_tensors with non-zero components providing fewer arguments")
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: PartExtracter does not know how to handle list_tensors with non-zero components providing fewer arguments

Process finished with exit code 1

Using Implict Euler or Explicit Euler works perfectly, I am not sure why the combination mix of the two doesn’t work?

Here is my code which works when timeIntScheme ='IE' or 'EE' for Implict or Explicit Euler respectively, but does not work when timeIntScheme ='SE' for Symplectic Euler

from fenics import *
import numpy as np

if __name__ == '__main__':

    timeIntScheme = 'IE' # WORKS
    # timeIntScheme = 'EE' # WORKS
    # timeIntScheme = 'SE' # DOESNT WORK

    # domain params
    nx = 1
    ny = 1
    xLength = 1.0
    yLength = 0.25

    # final time and steps
    tFinal = 20
    numSteps = 2000

    # ------------------------------# Create Mesh #-------------------------------#
    mesh = RectangleMesh(Point(xLength, 0.0), Point(0.0, yLength), nx, ny)

    # time and time step
    dt_value = tFinal/numSteps
    dt = Constant(dt_value)

    # long as the boundary between the two models are the same
    Bl_em = 1
    R1_em = 0.0
    R2_em = 0.0
    K_em = 1
    m_em = 2
    L_em = 0.1

    J_em = np.array([[0, Bl_em, 0], [-Bl_em, 0, -1], [0, 1, 0]])

    R_em = np.array([[R1_em, 0, 0], [0, R2_em, 0], [0, 0, 0]])

    A_conv_em = np.array([[1/L_em, 0, 0], [0, 1/m_em, 0], [0, 0, K_em]])

    A_em = Constant((J_em - R_em).dot(A_conv_em))

    # ------------------------------# Create Function Spaces and functions#-------------------------------#

    # create function spaces for mesh and for Electromechanical ODE variables
    odeSpace = FiniteElement('R', triangle, 0) #order 0, 1 variable
    element = MixedElement([odeSpace, odeSpace, odeSpace])
    U = FunctionSpace(mesh, element)

    # Define trial and test functions
    xode_0, xode_1, xode_2  = TrialFunctions(U)
    xode = as_vector([xode_0, xode_1, xode_2])
    v_xode_0, v_xode_1, v_xode_2 = TestFunctions(U)
    v_xode = as_vector([v_xode_0, v_xode_1, v_xode_2])

    # Define functions for solutions at previous and current time steps
    u_n = Function(U)
    xode_0_n, xode_1_n, xode_2_n = split(u_n)
    xode_n = as_vector([xode_0_n, xode_1_n, xode_2_n])

    u_ = Function(U)
    xode_0_, xode_1_, xode_2_ = split(u_)

    # -------------------------------# Boundary Conditions #---------------------------------#

    # Left edge boundary condition for marking
    class LeftMarker(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[0], 0.0)

    # Initialise mesh function for neumann boundary. Facets are dim()-1, initialise subdomains to index 0
    boundaryDomains = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
    boundaryDomains.set_all(0)

    leftMark= LeftMarker()
    leftMark.mark(boundaryDomains, 0)

    # redefine ds so that ds[0] is the dirichlet boundaries and ds[1] are the neumann boundaries
    ds = Measure('ds', domain=mesh, subdomain_data=boundaryDomains)

    # Forced input
    uInput = Expression('sin(3.14159*t)', degree=2, t=0)

    # -------------------------------# Problem Definition #---------------------------------#

    if timeIntScheme == 'SE':
        # THIS DOESNT WORK
        #implement Symplectic Euler integration
        xode_int = as_vector([xode[0], xode[1], xode_n[2]])
    elif timeIntScheme == 'EE':
        # THIS WORKS
        #implement Explicit Euler integration
        xode_int = as_vector([xode_n[0], xode_n[1], xode_n[2]])
    elif timeIntScheme == 'IE':
        # THIS WORKS
        #implement Implicit Euler integration
        xode_int = as_vector([xode[0], xode[1], xode[2]])

    # Create residual function
    F = (-dot(v_xode, xode - xode_n)/dt + dot(v_xode, A_em*xode_int) + \
             v_xode[0]*uInput )*ds(0)

    a = lhs(F)
    L = rhs(F)

    # Assemble matrices
    A = assemble(a)
    B = assemble(L)

    # Create vector for output displacement
    disp_vec = [0]

    # -------------------------------# Solve Loop #---------------------------------#

    t = 0
    for n in range(numSteps):
        set_log_level(LogLevel.ERROR)
        t += dt_value
        uInput.t = t

        B = assemble(L)

        #solve for this time step
        solve(A, u_.vector(), B)

        #output displacement
        disp = assemble(xode_2_*ds(0))/yLength # divide by yLength because 
#I integrate along the boundary length
        disp_vec.append(disp)

        # Update previous solution
        u_n.assign(u_)

Does anyone know why I am getting this error?

p.s This ODE will eventually be coupled with the boundary of a PDE. If anyone knows a better way to implement an ODE in Fenics I would be delighted to know!

Many thanks!

[Finbar_Argus]

This problem is fixed if I write out each term of dot(v_xode, A_em*xode_int) separately.

# This doesn't work
F = (-dot(v_xode, xode - xode_n)/dt + dot(v_xode, A_em*xode_int) + \
             v_xode[0]*uInput )*ds(0)

# This also doesn't work
F = (-dot(v_xode, xode - xode_n)/dt + dot(v_xode, dot(A_em, xode_int)) + \
             v_xode[0]*uInput )*ds(0)

# This works
F =(-dot(v_xode, xode - xode_n)/dt + \
      v_xode[0]*(A_em[0, 0]*xode_int[0] + A_em[0, 1]*xode_int[1] + A_em[0, 2]*xode_int[2]) + \
      v_xode[1]*(A_em[1, 0]*xode_int[0] + A_em[1, 1]*xode_int[1] + A_em[1, 2]*xode_int[2]) + \
      v_xode[2]*(A_em[2, 0]*xode_int[0] + A_em[2, 1]*xode_int[1] + A_em[2, 2]*xode_int[2]) +\
      v_xode[0]*uInput )*ds(0)

I have written out what dot( should do and there is no error, so I don’t know why using dot( causes an error?