Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Dx not working with 1D mixed function space #100

Open
j-bowhay opened this issue Nov 18, 2023 · 4 comments
Open

Dx not working with 1D mixed function space #100

j-bowhay opened this issue Nov 18, 2023 · 4 comments

Comments

@j-bowhay
Copy link

I am trying to solve the 1D Brusselator Reaction-diffusion equation:
image

However I am having issues trying to take the derivative of my trial function. Here is my code:

from sympy import symbols
from sympy.functions import cos, sin
from mpi4py_fft import generate_xdmf
from shenfun import inner, div, grad, TestFunction, TrialFunction, Function, \
    ETDRK4, FunctionSpace, Array, \
    comm, get_simplified_tpmatrices, Dx, MixedFunctionSpace

# Use sympy to set up initial condition
x = symbols("x", real=True)
u_1_0 = cos(x)
u_2_0 = sin(x)

# Initial conditions
u0 = 20*cos(x)
v0 = 20*sin(x)

# Parameters
a = 2
b = 6
D = 10

# Size of discretization
T1 = FunctionSpace(100, 'F', dtype='D', domain=(0, 20))
T2 = FunctionSpace(100, 'F', dtype='D', domain=(0, 20))
T = MixedFunctionSpace([T1, T2])
u = TrialFunction(T)
v = TestFunction(T)
padding_factor = 1.5

# Declare solution arrays and work arrays
U = Array(T, buffer=(u_1_0, u_2_0))
U_hat = Function(T)
w0 = Function(T)       # Work array spectral space
w1 = Array(T)         # Work array physical space

#initialize
U_hat = U.forward(U_hat)

def LinearRHS(self, u, dt):
    return Dx(u, 1, 2)

def NonlinearRHS(self, uv, uv_hat, rhs):
    rhs.fill(0)
    UVp = uv_hat.backward(padding_factor=padding_factor) # 3/2-rule dealiasing for nonlinear term
    w1[0] = a + UVp[0]**2*UVp[1] - (b+1)*UVp[0]
    w1[1] = b*UVp[0] - UVp[0]**2*UVp[1]
    rhs += w1.forward(w0)
    return rhs

if __name__ == '__main__':
    dt = 0.1
    end_time = 10
    integrator = ETDRK4(T, L=LinearRHS, N=NonlinearRHS)
    integrator.setup(dt)
    UV_hat = integrator.solve(U, U_hat, dt, (0, end_time))

However I get the following error:

Traceback (most recent call last):
  File "/home/jakeb/development/Further-Math-Bio-Special-Topic/code/1dspect.py", line 54, in <module>
    integrator.setup(dt)
  File "/home/jakeb/miniconda3/envs/shenfun/lib/python3.12/site-packages/shenfun/utilities/integrators.py", line 438, in setup
    L = self.LinearRHS(u, **self.params)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/jakeb/development/Further-Math-Bio-Special-Topic/code/1dspect.py", line 40, in LinearRHS
    return Dx(u, 1, 2)
           ^^^^^^^^^^^
  File "/home/jakeb/miniconda3/envs/shenfun/lib/python3.12/site-packages/shenfun/forms/operators.py", line 204, in Dx
    test = Dx(test, x, 1)
           ^^^^^^^^^^^^^^
  File "/home/jakeb/miniconda3/envs/shenfun/lib/python3.12/site-packages/shenfun/forms/operators.py", line 227, in Dx
    sc0 = sp.simplify(sp.diff(sc[i][j], psi[x], 1), measure=coors._measure)
                                        ~~~^^^
IndexError: tuple index out of range

This seems to be a bug(?) in the way Dx works as my trial function has two components however the coordinates of my function space are only 1d so psi only has one component so there is an index error in the loop here:

for i in range(dtest.num_components()):
for j in range(num_terms[i]):
sc0 = sp.simplify(sp.diff(sc[i][j], psi[x], 1), measure=coors._measure)

@mikaem
Copy link
Member

mikaem commented Nov 20, 2023

Hi
Unfortunately, the integrators do not work for a MixedFunctionSpace in 1D. However, since it is 1D it is not very difficult to implement a timestepper on your own using a BlockMatrixSolver. See below.

The MixedFunctionSpace has not really been kept up to date and should probably be used with care. I noticed that it did not even have a get_dealiased function and I have now added that to the master branch in order for this to work. I need to fix it even more, but with dealiasing in place (you need the master branch) I think you can implement your equations with something like:


from sympy import symbols, pi
from sympy.functions import cos, sin
from shenfun import inner, div, grad, TestFunction, TrialFunction, Function, \
    ETDRK4, FunctionSpace, Array, \
    comm, get_simplified_tpmatrices, Dx, MixedFunctionSpace, la

# Use sympy to set up initial condition
x = symbols("x", real=True)
u_1_0 = cos(x*pi)
u_2_0 = sin(x*pi)

# Initial conditions
u0 = 20*cos(x*pi)
v0 = 20*sin(x*pi)

# Parameters
a = 2
b = 6
D = 10

# Size of discretization
T1 = FunctionSpace(100, 'F', dtype='D', domain=(0, 20))
T2 = FunctionSpace(100, 'F', dtype='D', domain=(0, 20))
T = MixedFunctionSpace([T1, T2])
u = TrialFunction(T)
v = TestFunction(T)
padding_factor = 1.5

# Declare solution arrays and work arrays
UV = Array(T, buffer=(u_1_0, u_2_0))
UV_hat = Function(T)
w0 = Function(T)       # Work array spectral space
Tp = T.get_dealiased(padding_factor=1.5)
w1 = Array(Tp)

#initialize
UV_hat = UV.forward(UV_hat)

u0, u1 = u
v0, v1 = v
dt = 0.01
A = inner(v0, u0 - dt*Dx(u0, 0, 2)) + inner(v1, u1-dt*Dx(u1, 0, 2))
sol = la.BlockMatrixSolver(A)

def NonlinearRHS(uv_hat, rhs):
    rhs.fill(0)
    UVp = uv_hat.backward(padding_factor=padding_factor) # 3/2-rule dealiasing for nonlinear term
    w1[0] = a + UVp[0]**2*UVp[1] - (b+1)*UVp[0]
    w1[1] = b*UVp[0] - UVp[0]**2*UVp[1]
    rhs += w1.forward(rhs)
    return rhs

if __name__ == '__main__':
    rhs = Function(T)
    Nt = 4
    for n in range(Nt):
        rhs[:] = NonlinearRHS(UV_hat, rhs)
        UV_hat = sol(UV_hat + dt*rhs, u=UV_hat)

@j-bowhay
Copy link
Author

Thanks for the quick reply @mikaem! Would you kindly be able to explain the steps

A = inner(v0, u0 - dt*Dx(u0, 0, 2)) + inner(v1, u1-dt*Dx(u1, 0, 2))
sol = la.BlockMatrixSolver(A)

in the implementation, is this some form of backwards euler?

@mikaem
Copy link
Member

mikaem commented Nov 20, 2023

Yes, this is supposed to be Backwards Euler. But please note I've not tested it. I actually just realised that the last line should be

UV_hat = sol(dt*rhs, u=UV_hat)

and not

UV_hat = sol(UV_hat + dt*rhs, u=UV_hat)

because the UV_hat is already in A.

@mikaem
Copy link
Member

mikaem commented Nov 20, 2023

So sorry, got confused. Please disregard the last comment. It was correct in the first place.

Please note that you can also put parts of the (now) nonlinear term in the linear A.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants