[deal.II] Computing DFT

[Praveen C]

Dear all

I want to compute the DFT of a PDE solution. I have seen some Fourier Series in 
deal.II 

https://www.dealii.org/developer/doxygen/deal.II/classFESeries_1_1Fourier.html

and it says it is used to compute the DFT on the reference cell. I have seen 
examples of its use to estimate local  smoothness of solutions. 

Can these functions be used to compute the DFT on the whole domain, e.g., to 
compute things like energy spectrum of a turbulent velocity field ?

If this is not possible with these functions, do you have any suggestions how 
we can compute such spectra ?

Thanks
praveen

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/5C834555-DA13-4B32-AAC5-3E14FE5180FA%40gmail.com.

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

[Sylvain Mathonnière]

 Thank you for the answer. I have been looking at a way of doing what you 
suggested and I found the function* FEInterfaceValues 
::interface_dof_to_dof_indices()*
 
which seems to be doing what I want. 
I provide it with the interface index "i" and it returns a pair of local 
indices corresponding to the active and neighbouring cell. If then I query 
the first element  *dof_index[0]*, I obtain the local index of the current 
cell, so I though of using it like this :

*//get the local indices of the current cell and neighbouring 
cell*
*std::array dof_index = 
fe_iv.interface_dof_to_dof_indices(i);*
 
*   //use local index to obtain component of shape function.*
*   unsigned int component_i = 
fe_RTE.system_to_component_index(dof_index[0]).first;*

Then with the debugger and printing some values I realised that 
dof_index[0] is kind of always below 160 as expected, but on a couple of 
occasion it goes to 22116; 89130 or 4294967295 and triggers the error of 
system_to_component_index(). Those values looks like global indices whereas 
interface_dof_to_dof_indices(i) should returns local_indices. I must be 
doing something horribly wrong but I could not find a tutorial that uses 
interface_dof_to_dof_indices(). Am I confused with how to use those 
functions ?

Best,

Sylvain

El jueves, 22 de julio de 2021 a la(s) 21:53:19 UTC+2, Wolfgang Bangerth 
escribió:

>
> > Here *n_dofs = 320* and so *system_to_component_index()* does not like 
> to have 
> > an argument i>=160 and throws an error.
> > How come n_dofs=320 in this case ? I was expecting 2 interfaces * 2 
> nodes * 40 
> > components = 160 dofs, unless I am considering all the interfaces of the 
> cell 
> > at onces in which case I have 2 interfaces * 4 nodes * 40 components = 
> 320.
> > 
> > *My main question is therefore:*
> > How can I access the ith component of my shape function in face_worker 
> like in 
> > the cell_worker if I cannot use system_to_component_index() ?
>
> You are using FEInterfaceValues, which considers the shape functions that 
> live 
> on the interface. This set of shape functions is the union of the shape 
> functions living on the two adjacent cells. Because you have a 
> discontinuous 
> element, this set has size 2 * the number of dofs per cell, which gives 
> you 
> exactly the 2*160=320 you observe.
>
> Since you can't use FiniteElement::system_to_component_index on your index 
> 'i', you need to find a different way. One approach is to compute a 
> mapping 
> from this index 'i' to a pair '(here_or_there, index within cell)' and 
> then 
> call system_to_component_index() on the second half of the pair.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/2dcbdfa0-8a7f-4520-8614-558df5cf2287n%40googlegroups.com.

Re: [deal.II] Computing DFT

[Wolfgang Bangerth]

On 7/23/21 3:00 AM, Praveen C wrote:

I want to compute the DFT of a PDE solution. I have seen some Fourier Series in 
deal.II

https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fdeveloper%2Fdoxygen%2Fdeal.II%2FclassFESeries_1_1Fourier.html&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cc25d3f80dfac4442987708d94db85ce4%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637626276513834549%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=b3URzOQcYmPL41y86tT7g0eJcfW7oNu2rj%2Bcnf4TVdQ%3D&reserved=0

and it says it is used to compute the DFT on the reference cell. I have seen 
examples of its use to estimate local  smoothness of solutions.

Can these functions be used to compute the DFT on the whole domain, e.g., to 
compute things like energy spectrum of a turbulent velocity field ?

If this is not possible with these functions, do you have any suggestions how 
we can compute such spectra ?


Those functions are really just meant for the solution on a single cell. You 
will want to use one of the FFT packages out there to compute whole-domain 
spectra, say FFTx or FFTW.


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/9f750f68-0a54-d870-aad4-06111e8b0d95%40colostate.edu.

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

[Wolfgang Bangerth]

On 7/23/21 8:59 AM, Sylvain Mathonnière wrote:
Thank you for the answer. I have been looking at a way of doing what you 
suggested and I found the function*FEInterfaceValues 
::interface_dof_to_dof_indices()* 
which seems to be doing what I want.


Ah yes, I had forgotten about this function!


I provide it with the interface index "i" and it returns a pair of local 
indices corresponding to the active and neighbouring cell. If then I query the 
first element *dof_index[0]*, I obtain the local index of the current cell, so 
I though of using it like this :


*    //get the local indices of the current cell and neighbouring cell*
**
*    std::array dof_index = 
fe_iv.interface_dof_to_dof_indices(i);*

**
**
**
*   //use local index to obtain component of shape function.*
**
*   unsigned int component_i = 
fe_RTE.system_to_component_index(dof_index[0]).first;*


Then with the debugger and printing some values I realised that dof_index[0] 
is kind of always below 160 as expected, but on a couple of occasion it goes 
to 22116; 89130 or 4294967295 and triggers the error of 
system_to_component_index(). Those values looks like global indices whereas 
interface_dof_to_dof_indices(i) should returns local_indices. I must be doing 
something horribly wrong but I could not find a tutorial that uses 
interface_dof_to_dof_indices(). Am I confused with how to use those functions ?


No, I think you're right that that is how the function is supposed to work.

4294967295 is numbers::invalid_dof_index, which is the case discussed in the 
documentation where a DoF lives only one one side of the interface. The other 
values (22116; 89130) shouldn't happen. If that's what you get, then there is 
a bug somewhere.


Do you think you could construct a small testcase that illustrates what you 
have?

Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/c8c07ffd-3088-ba69-b911-5f775135fbf9%40colostate.edu.

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

[Sylvain Mathonnière]

 Oh, I think I know what those values appears and it is my bad ...
It is probably not 22116 but 22 and 116 and it is probably not 89130 but 
rather 89 and 130. I used something like *std::cout << dof_index[0] << 
std::endl* to get those values and sometimes they just get attached to each 
other like if the *std::endl* was skipped (not sure why, maybe because I 
have several threads ?). I did not realise it before reading your message 
but it actually makes sense because I was not getting an error for those 
values and my Assert was not triggered. I should have seen it, sorry.

As for the 4294967295 occurence, is it because I am dealing with an 
interface which is a boundary as well ? Isn't this case handled by the 
*boundary_worker* and not the *face_worker* (using step-12 notation) ?

Best,

Sylvain

El viernes, 23 de julio de 2021 a la(s) 17:17:17 UTC+2, Wolfgang Bangerth 
escribió:

> On 7/23/21 8:59 AM, Sylvain Mathonnière wrote:
> > Thank you for the answer. I have been looking at a way of doing what you 
> > suggested and I found the function*FEInterfaceValues 
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdealii.org%2Fdeveloper%2Fdoxygen%2Fdeal.II%2FclassFEInterfaceValues.html&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cbe067ec5c8ab4601e70308d94dea8304%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637626491881337201%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=hrnX9t7udM50eZvFVTiCBunh8ShGkow1iyI4AZq7udw%3D&reserved=0>::interface_dof_to_dof_indices()*
>  
>
> > which seems to be doing what I want.
>
> Ah yes, I had forgotten about this function!
>
>
> > I provide it with the interface index "i" and it returns a pair of local 
> > indices corresponding to the active and neighbouring cell. If then I 
> query the 
> > first element *dof_index[0]*, I obtain the local index of the current 
> cell, so 
> > I though of using it like this :
> > 
> > *//get the local indices of the current cell and 
> neighbouring cell*
> > **
> > *std::array dof_index = 
> > fe_iv.interface_dof_to_dof_indices(i);*
> > **
> > **
> > **
> > *   //use local index to obtain component of shape function.*
> > **
> > *   unsigned int component_i = 
> > fe_RTE.system_to_component_index(dof_index[0]).first;*
> > 
> > Then with the debugger and printing some values I realised that 
> dof_index[0] 
> > is kind of always below 160 as expected, but on a couple of occasion it 
> goes 
> > to 22116; 89130 or 4294967295 and triggers the error of 
> > system_to_component_index(). Those values looks like global indices 
> whereas 
> > interface_dof_to_dof_indices(i) should returns local_indices. I must be 
> doing 
> > something horribly wrong but I could not find a tutorial that uses 
> > interface_dof_to_dof_indices(). Am I confused with how to use those 
> functions ?
>
> No, I think you're right that that is how the function is supposed to work.
>
> 4294967295 is numbers::invalid_dof_index, which is the case discussed in 
> the 
> documentation where a DoF lives only one one side of the interface. The 
> other 
> values (22116; 89130) shouldn't happen. If that's what you get, then there 
> is 
> a bug somewhere.
>
> Do you think you could construct a small testcase that illustrates what 
> you have?
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/55622014-8692-448a-92ab-818c134c920an%40googlegroups.com.

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

[Corbin Foucart]

Sylvain,

Oh, I think I know what those values appears and it is my bad ...
> It is probably not 22116 but 22 and 116 and it is probably not 89130 but 
> rather 89 and 130. I used something like *std::cout << dof_index[0] << 
> std::endl* to get those values and sometimes they just get attached to 
> each other like if the *std::endl* was skipped (not sure why, maybe 
> because I have several threads ?). I did not realise it before reading your 
> message but it actually makes sense because I was not getting an error for 
> those values and my Assert was not triggered. I should have seen it, sorry.
>

 
I ran into this exact issue a few weeks ago. For printing statements while 
debugging in assembly loops, check the signatures of functions like 
MeshWorker::run and MeshWorker::mesh_loop. You should be able to pass 
parameters that force non-multithreaded execution, and your print 
statements will make sense. I found this helpful for debugging.

Best,
Corbin

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/7ec902ce-e1fa-4baa-a669-4b2b288479dfn%40googlegroups.com.

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

[Wolfgang Bangerth]

On 7/23/21 9:41 AM, Sylvain Mathonnière wrote:


As for the 4294967295 occurence, is it because I am dealing with an interface 
which is a boundary as well ? Isn't this case handled by the *boundary_worker* 
and not the *face_worker* (using step-12 notation) ?


No, it's because you are using DG elements. I tried to clarify this here:
  https://github.com/dealii/dealii/pull/12595

Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/7b60ecfa-9460-7339-91b3-2a4d1811ddda%40colostate.edu.

Re: [deal.II] Computing DFT

[Praveen C]

Thank you.

Related to this question.

Can I save a parallel distributed solution using serialize and read it back as 
a serial triangulation/solution ? This may  help me to sample the solution on a 
Cartesian grid and do FFT.

Thanks
praveen

> On 23-Jul-2021, at 8:30 PM, Wolfgang Bangerth  wrote:
> 
> Those functions are really just meant for the solution on a single cell. You 
> will want to use one of the FFT packages out there to compute whole-domain 
> spectra, say FFTx or FFTW.
> 
> Best
> W.

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/E623378D-BA62-4167-B15F-02BA90288E68%40gmail.com.