[deal.II] Re: DG explicit time integration for linear advection equation with MeshWorker (suggestions)

2019-09-18 Thread vachan potluri 


step-33 does compute interior face integrals twice.

One way I handle this is to attach a cell user index

 Thanks Praveen, this is similar to owner/neighbour concept of OpenFOAM :).


Doug,
Many thanks for the detailed explanations and your code :).

-- 
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/9b315e3b-c862-4e8f-8f78-c572cedc1f82%40googlegroups.com.

Re: [deal.II] Re: DG explicit time integration for linear advection equation with MeshWorker (suggestions)

2019-09-17 Thread Praveen C 


> On 18-Sep-2019, at 11:40 AM, Doug  wrote:
> 
> Praveen,
> 
> I did inspire my loops out of your dflo ;). Note that if you use a 
> distributed grid, the user_index() will not be consistent across MPI 
> boundaries. Therefore, you might end up with doubly-counted faces. Here is a 
> [commit](https://github.com/dougshidong/PHiLiP/commit/d0d7e6bd663b0cef032c7f953d2b2dddce4950a7#diff-d2da12f45bddf4e965acafc118ff8366
>  
> )
>  to your if with an if-statement that works for refine/distributed grids.

Thanks a lot for pointing this out. I suppose I never stumbled across this 
issue since I used meshworker in my parallel code.

Best regards
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/1582FBD1-6C12-4E28-876A-995BCBB0F227%40gmail.com.

[deal.II] Re: DG explicit time integration for linear advection equation with MeshWorker (suggestions)

2019-09-17 Thread Doug 
Hello Vachan,

Quick note. I'm only suggesting looping yourself because I am not totally 
familiar with MeshWorker. Looping myself exposed the mechanics a bit more 
for me, which I wanted.

I was pointing to step-33 for the loop, not the operators. I don't remember 
seeing a deal.II example implementing deal.II in operator form. 

So, I suggested looping though the cells to compute the operators that are 
not common to all matrices. This is usually the mass matrix if your cells 
are not linear (aka cartesian mesh). The other operators can be made the 
same for all the cells, except for some cofactor Jacobian matrix to 
transfer derivatives and normals from reference space to physical space. If 
you only ever plan on using a cartesian mesh, then all your local matrices 
will be the same. There is no point in assembling a global matrix that will 
operate on your solution. If I were you, I'd take a look at chapter 3 & 6 
of Hesthaven's book where he builds the same local operators that will be 
used on each cell. Therefore, you wouldn't have an array of differentiation 
matrices. Just one small differentiation matrix that applies to all cells. 
Feel free to look at how I do the [mass matrix](
https://github.com/dougshidong/PHiLiP/blob/master/src/dg/dg.cpp#L752)

Step-33 is CG by the way, therefore, two cells at the same level (no 
hanging nodes) don't need to be added to the residual. Therefore, they 
don't show that case. Feel free to look at my [loops](
https://github.com/dougshidong/PHiLiP/blob/master/src/dg/dg.cpp#L246), 
where I loop over all the cells, and all the faces. For internal faces, 
only one cell will be responsible to evaluate the flux and add it to both 
sides.

If someone wants to chime in with a better way to loop over all the faces 
just once, I would love to know about it too.

You can also keep on eye out on the code I linked. It currently does not 
pre-compute the operators, but it's only a matter of moving the loops out 
to pre-processing. 

Praveen,

I did inspire my loops out of your dflo ;). Note that if you use a 
distributed grid, the user_index() will not be consistent across MPI 
boundaries. Therefore, you might end up with doubly-counted faces. Here is 
a [commit](
https://github.com/dougshidong/PHiLiP/commit/d0d7e6bd663b0cef032c7f953d2b2dddce4950a7#diff-d2da12f45bddf4e965acafc118ff8366)
 to 
your if with an if-statement that works for refine/distributed grids.

On Wednesday, September 18, 2019 at 12:33:43 AM UTC-4, vachan potluri wrote:
>
> Doug and Praveen,
>
> Thanks for your answers. I had a look at step-33. As far as I understand, 
> although the looping through cells is not through MeshWorker, the assembly 
> is still global! So, for a non-cartesian mesh, I think you are suggesting 
> using such a loop over all cells to calculate local matrices. Will these 
> cell-local matrices stored as an array of matrices in the conservation law 
> class?
>
> I now have one more question. In assemble_face_term function of step-33, 
> the normal numerical flux is calculated. This function is called for every 
> cell-neighbor pair from the assemble_system function. This means, if I am 
> not wrong, at every interior face quadrature point, the numerical flux is 
> calculated twice. This might be very costly. Is there a way to avoid this? 
> One can probably force only one cell to calculate numerical flux at a face 
> based on the face normal's orientation. But then how can we communicate the 
> calculated flux to the other cell. Or even better, is it possible to loop 
> over faces? We could then add contributions to both the cells sharing this 
> face without double computation.
>
> Thanks again!
>

-- 
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/6f1ae5d3-a150-457a-8a6d-fa393910b15e%40googlegroups.com.

Re: [deal.II] Re: DG explicit time integration for linear advection equation with MeshWorker (suggestions)

2019-09-17 Thread Praveen C 


> I now have one more question. In assemble_face_term function of step-33, the 
> normal numerical flux is calculated. This function is called for every 
> cell-neighbor pair from the assemble_system function. This means, if I am not 
> wrong, at every interior face quadrature point, the numerical flux is 
> calculated twice. This might be very costly. Is there a way to avoid this? 
> One can probably force only one cell to calculate numerical flux at a face 
> based on the face normal's orientation. But then how can we communicate the 
> calculated flux to the other cell. Or even better, is it possible to loop 
> over faces? We could then add contributions to both the cells sharing this 
> face without double computation.
step-33 does compute interior face integrals twice.

One way I handle this is to attach a cell user index 

// set cell number
unsigned int index=0;
for (typename Triangulation::active_cell_iterator
cell=triangulation.begin_active();
cell!=triangulation.end(); ++cell, ++index)
cell->set_user_index (index);

and calculate face integral only once

for (unsigned int f=0; f::faces_per_cell; ++f)
{

if(cell->face(f)->at_boundary)
{
// boundary flux integral
}
else if(cell->neighbor_is_coarser(f))
{
// assemble from finer side
}
else if(cell->user_index() < cell->neighbor(f)->user_index())
{
// compute face integral
}

// add residual to left cell
// subtract residual from right cell
}

You add/subtract face integral to both cells adjacent to the face.

If you have adapted grid, then you assemble from the finer side. 

Best
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/918D8227-43EF-4CB6-A66E-234A077B9E5E%40gmail.com.

[deal.II] Re: DG explicit time integration for linear advection equation with MeshWorker (suggestions)

2019-09-17 Thread vachan potluri 
Doug and Praveen,

Thanks for your answers. I had a look at step-33. As far as I understand, 
although the looping through cells is not through MeshWorker, the assembly 
is still global! So, for a non-cartesian mesh, I think you are suggesting 
using such a loop over all cells to calculate local matrices. Will these 
cell-local matrices stored as an array of matrices in the conservation law 
class?

I now have one more question. In assemble_face_term function of step-33, 
the normal numerical flux is calculated. This function is called for every 
cell-neighbor pair from the assemble_system function. This means, if I am 
not wrong, at every interior face quadrature point, the numerical flux is 
calculated twice. This might be very costly. Is there a way to avoid this? 
One can probably force only one cell to calculate numerical flux at a face 
based on the face normal's orientation. But then how can we communicate the 
calculated flux to the other cell. Or even better, is it possible to loop 
over faces? We could then add contributions to both the cells sharing this 
face without double computation.

Thanks again!

-- 
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/92ef6323-88de-4dc2-8d80-72a0cbd617ec%40googlegroups.com.

[deal.II] Re: DG explicit time integration for linear advection equation with MeshWorker (suggestions)

2019-09-17 Thread Doug 
Hello Vachan,

What you are describing is very similar to Hesthaven's framework of build 
matrix operators and apply them to your degrees of freedom. The main 
advantage of this operator form is to have a common mass (if linear 
elements), differentiation, flux, and boundary operators for all your 
elements since the operators are defined on the reference element. 
Therefore, using MeshWorker loop to go through every single element to 
build those matrices would result in building redundant local matrices and 
storing them in a global matrix. Can be very expensive memory wise.

Yes, you could use MeshWorker to compute only your normal fluxes. 

So, it's feasible, but I wouldn't use the MeshWorker loop to build those 
matrices . You are welcome to look at step-33 to see how you can loop 
through your elements, and then build your different mass matrices. Then 
use a similar loop to apply your, interpolation, differentiation, lifting 
matrices on the local solution to assemble your right-hand side.

Yes, if you build those operators, you will save a significant amount of 
time *for explicit time stepping* since you will basically have 
concatenated multiple operators together. Interpolation, differentiation, 
integration, etc. into a single matrix-vector product.

I think know what you are referring to with that sparse flux matrix. If 
it's a local matrix, you would want to store the lifting operator, which is 
the mass inverse times your sparse flux matrix, which is not necessarily a 
sparse matrix (unless mass is diagonal). Either way, I wouldn't worry too 
much about that operation taking a long time. Your other operations on the 
volume take much longer, and if you decide to go with nonlinear fluxes, the 
time it takes to apply this lifted flux vector is meaningless. deal.II has 
a SparseMatrix class, in case you want to use it, and you can provide the 
set the SparsityPattern if you want. 

If you were planning on building multiple global matrices, then yes, your 
approach makes sense using SparseMatrix; be careful about memory. If you 
only plan on doing linear advection, your problems should be simple enough 
that you wouldn't care about memory nor computational time. I haven't seen 
people build those global matrices, but it should be faster than assembling 
every loop. It's the old trade-off between memory and computation.

Doug

On Tuesday, September 17, 2019 at 3:03:52 AM UTC-4, vachan potluri wrote:
>
> Hello all,
>
>
> I am a beginner in dealii. I want to solve a linear, transient advection 
> equation explicitly in two dimensions using DG. The resulting discrete 
> equation will have a mass matrix as the system matrix and a sum of terms 
> which depend on previous solution (multiplied by mass, differentiation, 
> flux and boundary matrix) as the rhs.
>
> [image: linearAdvection2D.png]
>
> Instead of using MeshWorker::loop for every single time step, I think the 
> following approach would be better. I am using a ghost cell approach to 
> specify the boundary condition: the boundary condition can be specified by 
> an appropriately calculated normal numerical flux.
>
>
>1. Before the any time steps, use a MeshWorker::loop each for each of 
>the four matrices: mass, differentiation, flux and boundary
>2. During each update
>   1. Again use MeshWorker::loop, but this time only to calculate the 
>   normal numerical flux.
>   2. Use the normal numerical flux and the previous solution to 
>   obtain the RHS using appropriate matrix-vector products
>   3. Solve the system
>
> I have few question regarding this approach.
>
>1. Is it feasible
>2. Can it give significant improvement in performance over the case 
>when assembling is done for every time step
>3. (Assuming answers to above questions are positive) For higher 
>orders, the flux and boundary matrices will be very sparse. The normal 
>numerical flux (which will be a vector) will also be sparse. Can the 
>matrix-vector product involving these combinations be optimised by using 
>appropriate sparsity pattern? Can a sparsity pattern be specified for a 
>vector too?
>
>

-- 
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/39f7baf9-e7c7-4a4a-b50e-b54ab5af0d7f%40googlegroups.com.