Simple integration doing something funny

[francispoulin]

Hello Mikeal,

I am trying to learn how to integrate arrays to be able to compute energies and other such scalars.

Below is some simple code that I tried to integrate 1 to hopefully get one but I find that I actually get 0.712. Clearly this is not working how I thought. Could you clarify?

Also, I want to do this on a vector space and wonder if there is a way to use inner cleverly or simply write it out componet wise?

Thanks in advance,
Francis

    V1  = FunctionSpace(N[0], 'F', dtype='D', domain=(0, L[0]))
    tmp = Array(V1, val=1.0)
    l2_error = np.linalg.norm(tmp,ord=2)/N[0]

[mikaem]

@francispoulin I'm not sure what you're trying to do, but did you see this https://mikaem.github.io/shenfun-demos/content/surfaceintegration.html ?

A 2D vector can be integrated with

Inner((1, 1), a)

Where 'a' is an 'Array' of a 2D 'VectorTensorProductSpace'

[francispoulin]

Thanks for the link and I had looked at that but you seemed to show me a prettier version than what I was reading.
Be;pw is a minimal example of some things that work and some things that fail.

The first print statement follows the example to integrate over the x-direction and produces the correct result.

The second print statement integrates over a 2D space and yields the correct result.

The third print statement tries to integrate a specific function but fails becuase 'C2C' object has no attribute 'bases'

The fourth print statement, the same as the first, does no work if you comment out the 3rd print statement. This seems to suggest that V1 changes somehow after we define T. Is that true?

    from mpi4py import MPI
    from shenfun import *
    from numpy import pi
    import sys

    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()

    M = 7
    N = (2**M, 2**M)
    L = [20.,  20.]

    V1  = FunctionSpace(N[0], 'F', dtype='D', domain=(0, L[0]))
    V2  = FunctionSpace(N[1], 'F', dtype='d', domain=(0, L[1]))

    # Works
    print(inner(1., Array(V1, val=1.0)).real/L[0])

    T   = TensorProductSpace(comm, (V1, V2), **{'planner_effort': 'FFTW_MEASURE'})
    TV  = VectorTensorProductSpace(T)

    X  = T.local_mesh(True)

    Ub =  Array(T)
    Ub =  1./pow(np.cosh(X[1] - L[1]/2),2)

    # Works
    print(inner(1., Array(T, val=1.0)).real/(L[0]*L[1]))

    # Fails
    print(inner(Ub, Array(T, val=1.0)).real/(L[0]*L[1]))

    # If this line is uncommented, the next line works
    #V1  = FunctionSpace(N[0], 'F', dtype='D', domain=(0, L[0]))
    print(inner(1., Array(V1, val=1.0)).real/L[0])

[francispoulin]

Just to try and share what I have learned, and what I don't understand, here is another sample code that does a few things that I don't quite understand.

    from mpi4py import MPI
    from shenfun import *
    from numpy import pi

    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()

    N = (8,8)
    L = (20., 20.)

    V1  = FunctionSpace(N[0], 'F', dtype='D', domain=(0, L[0]))
    V2  = FunctionSpace(N[1], 'F', dtype='d', domain=(0, L[1]))

    # Test 1D integrals
    print('Test 1:  Success, both are 1.0')
    print(inner(1, Array(V1, val=1)).real/L[0])
    print(inner(1, Array(V2, val=1)).real/L[1])
    print(' ')

    T   = TensorProductSpace(comm, (V1, V2), **{'planner_effort': 'FFTW_MEASURE'})
    TV  = VectorTensorProductSpace(T)

    # Test 2D integrals
    print('Test 2: Question, why are the results not the same? 2.0 vs 1.96875')
    print(inner( (1, 1), Array(TV, val=1))/(L[0]*L[1]))
    print(inner( (1, 1), Array(TV, val=(1,1)))/(L[0]*L[1]))
    print(' ')

    X  = T.local_mesh(True)

    Uj    = 1.0
    Ub    = Array(TV)
    Ub[0] = 1. + 0*X[1]
    Ub[1] = 0. + 0*X[1]

    print('Test 3: Question, why is this not the same as test 2? 9.17e+47')
    print(inner( (1, 1), Array(TV, val=(1,1)))/(L[0]*L[1]))
    print(' ')

    print('Test 4: Question, why does this syntax not work? AssertionError')
    print(inner( Ub, Array(TV, val=(1,1)))/(L[0]*L[1]))
    print(' ')

[mikaem]

Ok, a few things. This is not allowed:

print(inner(Ub, Array(T, val=1.0)).real/(L[0]*L[1]))

because an Array (Ub) is not a basis function. But this works:

print(inner(1, Ub*Array(T, val=1.0)).real/(L[0]*L[1]))

Here 1 is interpreted as a constant test function and the inner product becomes an unweighted integral of Ub*Array(T, val=1.0) over the domain.

#If this line is uncommented, the next line works
#V1 = FunctionSpace(N[0], 'F', dtype='D', domain=(0, L[0]))
print(inner(1., Array(V1, val=1.0)).real/L[0])

In the latest version of shenfun this works already. The reason is that in previous versions V1 was changed when it became a part of a multidimensional TensorProductSpace. In recent versions of shenfun a copy of V1 is made under the hood in the creation of the TensorProductSpace and V1 is unchanged. Thus the line print(inner(1., Array(V1, val=1.0)).real/L[0]) works as expected.

#Test 2D integrals
print('Test 2: Question, why are the results not the same? 2.0 vs 1.96875')
print(inner( (1, 1), Array(TV, val=1))/(L[0]*L[1]))
print(inner( (1, 1), Array(TV, val=(1,1)))/(L[0]*L[1]))

The reason is that I have not implemented this last syntax:

print(inner( (1, 1), Array(TV, val=(1,1)))/(L[0]*L[1]))

Look at the Array created. It contains gibberish. But the following syntax actually works

 print(inner( (1, 1), Array(TV, buffer=(1,1)))/(L[0]*L[1]))

I should definitely add your syntax as well:-)

print('Test 4: Question, why does this syntax not work? AssertionError')
print(inner( Ub, Array(TV, val=(1,1)))/(L[0]*L[1]))

Same as over. Works with inner( (1, 1), Ub*Array(TV, buffer=(1,1)).

[francispoulin]

Thank you Mikael, that was most helpful. It took me a bit of time to figure out (excuse previous post that I deleted) but I did get everything working from before. Great!

Unfortunately, I'm finding odd results with padded function spaces. I would expect, perhaps naively, even if I have a function space to be padded, it should yield the same results. They do not and I'm trying to figure this out. Any advice?

    # Padded Spaces
    Tp  =  T.get_dealiased((1.5, 1.5))
    TVp = TV.get_dealiased((1.5, 1.5))

    # Test 2D integrals, cannot use val yet for vector
    print('Test: Question, why do we get 2.0 and 4.5?')
    print(inner( (1, 1), Array(TV, buffer=(1,1)))/(L[0]*L[1]))
    print(inner( (1, 1), Array(TVp, buffer=(1,1)))/(L[0]*L[1]))
    print(' ')

    X, Xp  = T.local_mesh(True), Tp.local_mesh(True)

    Uj    = 1.0
    Ub, Ubp = Array(TV), Array(TVp)
    Ub[0], Ubp[0] = 1. + 0*X[1], 1. + 0*Xp[1]
    Ub[1], Ubp[1] = 0. + 0*X[1], 0. + 0*Xp[1]

    print('Test 4: Question, why do we get different values for padded functions? 1.0 vs 2.25')
    print(inner( 1, Ub[0]*Array(T, val=1))/(L[0]*L[1]))
    print(inner( (1,1), Ub*Array(TV, buffer=(1,1)))/(L[0]*L[1]))
    print(inner( 1, Ubp[0]*Array(Tp, val=1))/(L[0]*L[1]))
    print(inner( (1,1), Ubp*Array(TVp, buffer=(1,1)))/(L[0]*L[1]))
    print(' ')
 

[mikaem]

Good catch:-) I didn't consider these unweighted integrals with padding before. Frankly, I'm not 100% sure if it should be allowed, and any other space than Fourier would here throw an error. The padded arrays are there for computing non-linearities and nothing else really. If you want a finer resolution space there's the refine method. Also, you would not get higher accuracy from using padded arrays, it would simply be less efficient than a regular Array. So I think I prefer adding an error message if someone tries to integrate using a padded Array.

Minor reconsideration. If you compute some linearity on the padded array, then the accuracy of an integral over that array would be improved.

[francispoulin]

It sounds like maybe I shouldn't be considering this in the first place so maybe adding an error message is the simplest thing to do.

This came up because in my code I am defining my initial conditions on a padded mesh since I thought that's what you should be doing in order to have the aliasing work correctly. I will look at your sample codes and see what you do there and try and mimic that more closely.

[mikaem]

There's nothing wrong with initialising on a padded mesh, but the array and highest frequencies will be truncated moving to spectral space, so I don't think there's much of a benefit there.

[mikaem]

BTW - the fix for Fourier is a one-liner

Change [this line] (

shenfun/shenfun/fourier/bases.py

Line 109 in 919ae01

N = self.N

) with code

N = self.N

to

N = self.shape(False)

and everything works as expected:-) I'll check if there's any damage. If not, then I can use the fix.

[francispoulin]

Thanks for finding such an easy fix. I am happy to change that on my local machine but could you tell me where I would find shenfun installed by conda?

I think this has also pointed out that I should probably not be initializing my fields on padded cells, so it's good that I fix that. I am very glad to learn that now, just sorry that I made this silly mistake before.

[mikaem]

I just ran all tests with the implemented fix and it seems to be alright, so I'll keep it. You can find shenfun doing

import shenfun
print(shenfun.__file__)

[francispoulin]

Great! And I will try and not forget that information

I can confirm that it does work with this small and clever change. Thanks!