Comments (12)
@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'
from shenfun.
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])
from shenfun.
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(' ')
from shenfun.
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))
.
from shenfun.
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(' ')
from shenfun.
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.
from shenfun.
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.
from shenfun.
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.
from shenfun.
BTW - the fix for Fourier is a one-liner
Change [this line] (
shenfun/shenfun/fourier/bases.py
Line 109 in 919ae01
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.
from shenfun.
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.
from shenfun.
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__)
from shenfun.
Great! And I will try and not forget that information
I can confirm that it does work with this small and clever change. Thanks!
from shenfun.
Related Issues (20)
- Dx not working with 1D mixed function space HOT 4
- Issue on Stokes equations with doubly periodic boundary conditions. HOT 2
- Error when evaluating derivatives for nonstandard domain HOT 7
- Installation problem: error in "cimport fastgl_wrap" HOT 2
- Old shenfun script for 2D Poisson equation dosn't work with new version HOT 2
- 1d poisson HOT 1
- Bug in FunctionSpace with bc as dictionary HOT 4
- Solve DGL with non-constant (discrete) variable? HOT 5
- AttributeError: 'tuple' object has no attribute 'copy' HOT 1
- Install problem HOT 1
- problems in solving 2D Allen Cahn equation HOT 2
- Two questions in the demo MixedPoisson.py HOT 7
- How to compute the numerical solution on the general points HOT 4
- 1. Solver for multidimensional problem 2. Product with variable coefficients HOT 4
- Two questions about the Legendre-Galerkin method for the Cahn-Hilliard equation HOT 2
- How to compute the integral about the nonlinear term in the Navier-Stokes equation? HOT 2
- How to assemble the right-hand side vector? HOT 3
- Understanding the chebyshev solver in poisson2ND.py HOT 3
- Time dependent boundary conditions in segments for Rayleigh-Benard 2D HOT 2
- Implementing the one-dimensional Kuramoto-Sivashinsky system HOT 3
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from shenfun.