mirok / fenics-calc

Calculator over FEniCS functions

fenics-calc

This module provides a simple lazy evaluator of algebraic expressions over FEniCS functions. Expressions are defined in the UFL language. They can be evaluated provided that the function spaces of each Function type arguments are of the tensor product form V x V ... x V. The expression then results in a Function in some appropriate tensor product space with V. In general coefficients of the function for expression (op u v) are computed as (op U V) where U, V are coefficient vectors of the arguments. This approach means that the result is exact iff op is linear; otherwise there is an interpolation error.

from dolfin import *
from xcalc import Eval

mesh = UnitSquareMesh(5, 5)
T = TensorFunctionSpace(mesh, 'CG', 1)

A = interpolate(Expression((('x[0]', 'x[1]'),
                                    ('2*x[0]+x[1]', 'x[0]+3*x[1]')), degree=1), T)
expr = tr(sym(A) + skew(A))
me = Eval(expr)  # Function in VectorFunctionSpace(mesh, 'CG', 1)
                 # Exact due to linearity

Functions can be grouped into TempSeries (avoiding name collision with FEniCS's native TimeSeries). Same algebraic operations over these objects are supported as with normal functions. TempSeries are collapsed into functions by time-averaging operations such as mean

from dolfin import *
from xcalc.tempseries import PVDTempSeries
from xcalc.operators import Mean

# Let there be a time series stored in pvd of function in V
mesh = UnitSquareMesh(3, 3)
V = FunctionSpace(mesh, 'CG', 1)
series = PVDTempSeries('pvd_test.pvd', V)

mean = Mean(abs(series))  # A lazy node! 

What works and what does not

At the moment the interpreter supports most commonly used nodes in UFL except for

• differentiation nodes, e.g. grad (hence algebraic expressions)
• FEM specific nodes such as FacetNormal, jump, avg and so on
• Currently MPI support is missing.

Limitations

In order for Eval to be successfull operations on the dofs must translate nicely to numpy. Among other things this means that the spaces where the expression terminals are defined must be collapseable to scalars (in order to access components). This rules out vector-valued elements such as Raviart-Thomas

from dolfin import *
from xcalc import Eval

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'RT', 1)

x = Function(V)
y = Function(V)

 Eval(inner(x, y))  # Fails - what is the scalar space here?
 Eval(as_vector((x[0], y[1])))  # Fails - same reason
 Eval(2*x + y)  # Fails - handled internally as 2*as_vector((x[0], x[1])) + y
 Eval(dot(outer(x, y), x))  # Fails - numpy eval of the mat is array of length 2
                            #         numpy eval of the vec is single number
                            #         numpy op produces array of 2 where 1 is expected 

FEniCS compatibility

This package is CI tested against FEniCS packages for ubuntu 16.04 LTS

Dependencies and installation

In addition to FEniCS stack h5py is needed if one wants to use XDMFTempSeries. Dynamic mode decomposition of series relies on PyDMD module. After that, python setup.py install --user (or variants) is how to install this module.

fenics-calc's Issues

No module named 'interpreter'

@MiroK
Hi, Miroslav. Thanks for your work. I want to use this code to do some calculations, but I meet an error as :

Is there any requirement for the version of fenics ? I am using fenics 2018.1.0 (python3). I make a small modification to install fenics-xalc by python3.

Could you give me some suggestions, please ?

compatible version of PyDMD

@MiroK Hello, Miro. Do you know (or have a look the code you are using) the compatible version of PyDMD with your code ? I use following sample code, but it cannot run with latest version of PyDMD.

from xcalc.timeseries import PVDTempSeries
from xcalc.interpreter import Eval
from xcalc.dmd import dmd
from dolfin import *
import numpy as np

try:
    from pydmd import DMD
    # https://github.com/mathLab/PyDMD
except ImportError:  # Fallback to standard DMD if pydmd one is missing
    from xcalc.dmdbase import DMD

dmd_worker = DMD(svd_rank=-1)  # No truncation


what = 'pressure'
dir = 'ActiveControl'

if what == 'pressure':
    no_control = PVDTempSeries('%s/p_out.pvd' % dir, FiniteElement('Lagrange', triangle, 1))

elif what == 'velocity':
    # I chose to look here at velocity magniture, pod of velocity field shoudl work as well
    velocity = PVDTempSeries('%s/u_out.pvd' % dir, VectorElement('Lagrange', triangle, 1))
    
    no_control = Eval(sqrt(inner(velocity, velocity)))

# Descending according to energy
temp_slice = slice(200, 401, None)
# Look at time evolution of coefs for first 6  modes
which_modes = range(6)
spectrum, modes, dynamics = dmd(no_control.functions[temp_slice], dmd_worker, modal_analysis=which_modes)

print 'Reconstruction error', np.linalg.norm(dmd_worker.snapshots - dmd_worker.reconstructed_data.real)

rout = File('%s/dmd_%s_real_%s.pvd' % (dir, dir.lower(), what))
iout = File('%s/dmd_%s_imag_%s.pvd' % (dir, dir.lower(), what))

for index, f in enumerate(modes):
    # Plo real
    fr = f.real
    fr.rename('f', '0')
    rout << (fr, float(index))

    fi = f.imag
    fi.rename('f', '0')
    iout << (fi, float(index))

# Dump the energie
np.savetxt('%s/dmd_%s_%s_spectrum.txt' % (dir, dir.lower(), what), spectrum[which_modes])

import matplotlib
matplotlib.use('Agg')

import matplotlib.pyplot as plt

# Spectrum
fig = plt.figure()
ax = fig.gca()
plt.plot(spectrum.real, spectrum.imag, color='blue', marker='o', linestyle='none')

for i, (x, y) in enumerate(zip(spectrum.real, spectrum.imag)):
    plt.text(x, y, str(i))

c = plt.Circle((0., 0.), 1., color='green', fill=False, linestyle=':')  # Unit circle to see which modes blow up
ax.add_artist(c)
plt.axis('equal')
plt.savefig('%s/dmd_%s_%s_spectrum.png' % (dir, dir.lower(), what))

# Energy
fig = plt.figure()
ax = fig.gca()
plt.plot(np.abs(spectrum), color='blue', marker='o', linestyle='none')
plt.savefig('%s/dmd_%s_%s_energy.png' % (dir, dir.lower(), what))

# Evolution of modes
times = no_control.times[temp_slice]
indices = [611, 621, 631, 641, 651, 661]

fig, axarr = plt.subplots(6, 1)
for i, coef_i in enumerate(dynamics):
    axarr[i].plot(times, coef_i.real, label=str(i))
plt.savefig('%s/dmd_%s_%s_modes.png' % (dir, dir.lower(), what))

plt.show()