Define function with complicated expression and interpolate it

[cpraveen]

I want to interpolate a function which is defined as a sum of many terms.

Firstly this does not work and gives lot of errors in output.

It runs if I have a sum of 100 terms but not for 1000 terms.

Secondly, loop is going to be slow, how can I vectorize it ?

from firedrake import *

Ha = 10.0
Re = 10.0

def vexact(X):
    x, y = X[0], X[1]
    n = 1000
    K = Ha / (1 - 0.825 / sqrt(Ha) - 1.0 / Ha)
    dpdz = K / Re
    value = 0.0
    for i in range(n):
        alpha = (i + 0.5) * pi
        N = sqrt(Ha**2 + 4 * alpha**2)
        r1 = 0.5 * ( Ha + N )
        r2 = 0.5 * (-Ha + N )
        t1 = exp(-2 * r1)
        t2 = exp(-2 * r2)
        t3 = 1 - t1 * t2
        V2 = (1 - t2) * (exp(-r1 * (1 - y)) + exp(-r1 * (1 + y))) / (2 * t3)
        V3 = (1 - t1) * (exp(-r2 * (1 - y)) + exp(-r2 * (1 + y))) / (2 * t3)
        value += 2 * (-1)**i * (1 - V2 - V3) * cos(alpha * x) / alpha**3
    return dpdz * value

mesh = RectangleMesh(10, 10, 1.0, 1.0, originX=-1.0, originY=-1.0)

V = FunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)
vz = Function(V).interpolate(vexact(x))

[connorjward]

You are constructing an enormous symbolic expression, so I can see why this might crash and be very slow.
You might have more joy doing this numerically instead. mesh.coordinates.dat.data_ro will give you back a numpy array containing the mesh coordinates which you could try using instead of SpatialCoordinates. At the end you can put the resulting array into a Function by passing the val= keyword argument. To speed things up you might find numba helpful here.

[cpraveen]

Thank you, it is better now.

def vexact(x,y):
    n = 1000
    K = Ha / (1 - 0.825 / np.sqrt(Ha) - 1.0 / Ha)
    dpdz = K / Re
    i = np.arange(n)
    alpha = (i + 0.5) * np.pi
    N = np.sqrt(Ha**2 + 4 * alpha**2)
    r1 = 0.5 * ( Ha + N )
    r2 = 0.5 * (-Ha + N )
    t1 = np.exp(-2 * r1)
    t2 = np.exp(-2 * r2)
    t3 = 1 - t1 * t2
    V2 = (1 - t2) * (np.exp(-r1 * (1 - y)) + np.exp(-r1 * (1 + y))) / (2 * t3)
    V3 = (1 - t1) * (np.exp(-r2 * (1 - y)) + np.exp(-r2 * (1 + y))) / (2 * t3)
    value = 2 * (-1)**i * (1 - V2 - V3) * np.cos(alpha * x) / alpha**3
    return dpdz * np.sum(value)

mesh = RectangleMesh(100, 100, 1.0, 1.0, originX=-1.0, originY=-1.0)

V = FunctionSpace(mesh, "CG", 1)

# Interpolate velocity
xy = SpatialCoordinate(mesh)
x = Function(V).interpolate(xy[0])
y = Function(V).interpolate(xy[1])
x, y = x.dat.data_ro[:], y.dat.data_ro[:]

vzdata = np.zeros(len(x))
for i in range(len(x)):
    vzdata[i] = vexact(x[i],y[i])
vz = Function(V, val=vzdata, name="vz")