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

[vachan potluri]

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/7fd51309-f137-421f-91bc-8c2d1c6b7ff1%40googlegroups.com.

[deal.II] Problem warning in eclipse: invalid argument for the distance function

[Brian Zhang]

Hello, everyone,


I am pretty new to deal.ii.


I have studied the tutorial step-1, where it is recommended to learn a 
debugger ASAP. 


So I put it into the Eclipse IDE.


However, I found there is an error as shown in the bellow picture. It is 
related to that the argument of the "distance" function (distance(const 
dealii::Point &) should be const, while cell->vertex(v) is not 
const.


[image: Screenshot from 2019-09-17 15-36-45.png]


Even though the program can run successfully, I still feel uncomfortable 
with this error notification existing in the code.


So I made some modification as shown in the below. Now the error 
notification is gone.

[image: Screenshot from 2019-09-17 15-37-54.png]


This error notification looks very common in Eclipse IDE, for example 
another one about "set_manifold" function, as shown in the below screen 
shot. 

[image: Screenshot from 2019-09-17 15-33-18.png] 

The description of the problem has also been copied and pasted here: 
Invalid arguments ' Candidates are: 

void set_manifold(unsigned int, const ? &)
void set_manifold(unsigned int).



My question here is what is the routine for dealing with such problems when 
using Eclipse IDE for deal.ii? Ignore it or modify as what I have done?


Thank you very much for your time and help!


Brian Zhang

U of Waterloo

-- 
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/a162377f-4c51-4437-86f8-7bb1a2097e53%40googlegroups.com.

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

[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.

[deal.II] Instantiation MappingFEField with hp::DoFHandler

[Doug]

Hello again,

I am trying to use MappingFEField to represent my high-order mesh points 
(which will then be deformed through optimization). Similarly to this 
thread High order mesh from Gmsh 
.
 
However, I am using hp-finite element and I would expect different element 
geometric orders to exist for cells with different solution orders.

I see that MappingFEField is currently only instantiated with DoFHandler, 
and not hp::DoFHandler, so I tried to add it to the instantiation file 
"/fe/mapping_fe_field.inst.in". Unfortunately, it seems like the current 
MappingFEField implementation assumes that the finite element is not a 
FECollection, etc. and tries to directly operate on it. 

What would be the best way to have different element geometric orders, 
where I can control the metric degrees of freedom?

Doug

-- 
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/a97c4d3c-3557-40fa-b9fc-654bf3516551%40googlegroups.com.

Re: [deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[Doug]

>
> What I meant is: Can you show the compiler error message that illustrates 
> where the assertion is located, what the template arguments are, how it 
> came 
> that we called that function with these template arguments, etc? 
>
>
/home/ddong/Libraries/dealii/include/deal.II/base/numbers.h:583:3: note:  
 no known conversion for argument 1 from ‘const Sacado::Fad::DFad’ 
to ‘const std::complex&’
In file included from 
/home/ddong/Libraries/dealii/include/deal.II/base/cuda.h:21:0,
 from 
/home/ddong/Libraries/dealii/include/deal.II/base/memory_space.h:22,
 from 
/home/ddong/Libraries/dealii/include/deal.II/lac/la_parallel_vector.h:21,
 from 
/home/ddong/Libraries/dealii/source/lac/la_parallel_vector.cc:16:
/home/ddong/Libraries/dealii/include/deal.II/base/exceptions.h:1671:42: 
error: no matching function for call to 
‘std::complex::complex(const Sacado::Fad::DFad&)’
  dealii::ExcNumberNotFinite(std::complex(number)))

It is basically part of the same call to AssertIsFinite to generate the 
Exception. Note that it could hit this assertion much more directly 
since 
/home/ddong/Libraries/dealii/include/deal.II/lac/la_parallel_vector.templates.h 
itself calls AssertIsFinite. e.g. line 1450

In the end, it is basically trying to convert the Sacado type into a 
complex number, which is undefined, whenever it tries to perform some 
vector operations.

Doug

-- 
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/5fae8c77-687b-4127-b6b3-e62a05b24949%40googlegroups.com.

Re: [deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[Wolfgang Bangerth]

On 9/17/19 6:58 PM, Doug wrote:
> 
> /home/ddong/Libraries/dealii/include/deal.II/base/numbers.h:583:3: 
> note:   no known conversion for argument 1 from ‘const 
> Sacado::Fad::DFad’ to ‘const std::complex&’
> In file included from 
> /home/ddong/Libraries/dealii/include/deal.II/base/cuda.h:21:0,
>                   from 
> /home/ddong/Libraries/dealii/include/deal.II/base/memory_space.h:22,
>                   from 
> /home/ddong/Libraries/dealii/include/deal.II/lac/la_parallel_vector.h:21,
>                   from 
> /home/ddong/Libraries/dealii/source/lac/la_parallel_vector.cc:16:
> /home/ddong/Libraries/dealii/include/deal.II/base/exceptions.h:1671:42: 
> error: no matching function for call to 
> ‘std::complex::complex(const Sacado::Fad::DFad&)’
>            dealii::ExcNumberNotFinite(std::complex(number)))
> 
> It is basically part of the same call to AssertIsFinite to generate the 
> Exception. Note that it could hit this assertion much more directly 
> since 
> /home/ddong/Libraries/dealii/include/deal.II/lac/la_parallel_vector.templates.h
>  
> itself calls AssertIsFinite. e.g. line 1450

But I don't know how we get there. What does the rest of the error 
message look like? The place you show looks like this:

 template 
 void
 Vector::add(const Number a)
 {
   AssertIsFinite(a);

so 'a' is of type Sacado::Fad::DFad. It then needs to call

   numbers::is_finite (Sacado::Fad::DFad), which doesn't exist. 
It probably just tries to go through all of the overloads of 
numbers::is_finite() and wants to see whether it can convert the 
argument Sacado::Fad::DFad to the argument type of these 
overloads. The error message you show then would just explain why this 
one possibility (namely, numbers::is_finite(std::complex&)) 
does not work. But that doesn't mean that this is the right overload 
anyway -- I suspect that your compiler produces similar error messages 
above or below the one you show for all of the other overloads, right?

I *think* that the solution is to simply provide an overload for
   numbers::is_finite (const Sacado::Fad::DFad &x)
Can you try this? You could declare it in your own .cc file before you 
#include 

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/846b94a1-9c0c-babc-daa3-11653974365b%40colostate.edu.

Re: [deal.II] Cmake-GUI error

[Wolfgang Bangerth]

On 9/17/19 4:13 PM, GaryR wrote:
> I keep getting the following error when running cmake-gui.
> |
> 
> 
> CMakeErrorat 
> /home/gary/Deal.II/share/deal.II/macros/macro_deal_ii_setup_target.cmake:64(INCLUDE):
>    INCLUDE could notfind load file:
> 
> 
> /home/gary/Deal.II/lib/cmake/deal.II/deal.IITargets.cmake
> CallStack(most recent call first):
> /home/gary/Deal.II/share/deal.II/macros/macro_deal_ii_invoke_autopilot.cmake:55(DEAL_II_SETUP_TARGET)
> CMakeLists.txt:39(DEAL_II_INVOKE_AUTOPILOT)
> 
> |
> 
> I've tried mucking my way through the CMakeLists.txt. and the 
> setup_target_.cmake files. I can find the place where the error occurs 
> but am not familiar enough with the coding to make sense of the problem.
> 
> I am running Debian Buster with a KDE Desktop.
> AMD 4 core processor with 16GB of ram.
> Using Deal.II-9.1.1
> Using Cmake-GUI-3.13.2
> Using Eclipse 2019.6 (4.12.0)
> 
> I have Cmake-GUI set up to produce an Eclipse compatible build.

I don't think any of us uses the Cmake-GUI. We all just call cmake on 
the command line. Does that work?

It's also not entirely clear to me at which stage this happens to you. 
Is this the deal.II cmake you're trying to execute? Or in a tutorial 
program? Or in your own project?

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/ad0d5d16-4b76-ab0e-a3f5-b68e2364e67d%40colostate.edu.

Re: [deal.II] Problem warning in eclipse: invalid argument for the distance function

[Daniel Arndt]

Brian,

The warnings you are seeing come from a static analyzer in Eclipse and not
from the compiler. It seems like Eclipse does not have full access to the
source code
or doesn't understand it correctly. E.g., passing a non-const variable to a
function that doesn't modify its parameter, is never a problem.

It is certainly helpful to have a look if these warnings are sensible or
not. In this case, I would ignore them though. On the other hand, I would
listen to compiler warnings.

Best,
Daniel

-- 
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/CAOYDWbL5qjn6baPKRPBicfW36XGq-z7AL_%3DqCScNpzVkF6k8iQ%40mail.gmail.com.

Re: [deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[Doug]

>
> so 'a' is of type Sacado::Fad::DFad. It then needs to call 
>
>numbers::is_finite (Sacado::Fad::DFad), which doesn't exist. 
> It probably just tries to go through all of the overloads of 
> numbers::is_finite() and wants to see whether it can convert the 
> argument Sacado::Fad::DFad to the argument type of these 
> overloads. The error message you show then would just explain why this 
> one possibility (namely, numbers::is_finite(std::complex&)) 
> does not work. But that doesn't mean that this is the right overload 
> anyway -- I suspect that your compiler produces similar error messages 
> above or below the one you show for all of the other overloads, right? 
>
>
True, the error message gets long pretty quickly, but the undefined 
is_finite was another issue. Even if is_finite exists, the complex 
constructor is still an issue.
 

> I *think* that the solution is to simply provide an overload for 
>numbers::is_finite (const Sacado::Fad::DFad &x) 
> Can you try this? You could declare it in your own .cc file before you 
> #include  
>

 Attached is a tiny .cc file that simply instantiates the vector. Copy 
pasted here here convenience.

#include 
namespace dealii{
namespace numbers{
bool is_finite(const Sacado::Fad::DFad &x) {
¦   (void) x;
¦   return true;
}   
}}
#include 
#include 
int main (int /*argc*/, char * /*argv*/[]){
using namespace dealii;
dealii::LinearAlgebra::distributed::Vector vector_double;
using ADtype = Sacado::Fad::DFad;
dealii::LinearAlgebra::distributed::Vector vector_ad;
return 0;
}

I also provide an function for numbers::is_finite (const 
Sacado::Fad::DFad &x) to avoid the first set of errors. However, I 
still get the error below

In file included from 
/home/ddong/Libraries/dealii/install/include/deal.II/base/aligned_vector.h:22:0,
 from 
/home/ddong/Libraries/dealii/install/include/deal.II/lac/vector.h:22,
 from 
/home/ddong/Codes/PHiLiP/src/instantiate_vector_ad.cpp:13:
/home/ddong/Libraries/dealii/install/include/deal.II/lac/la_parallel_vector.templates.h:
 
In instantiation of ‘dealii::LinearAlgebra::distributed::Vector& dealii::LinearAlgebra::distributed::Vector::operator*=(Number) [with Number = Sacado::Fad::DFad; 
MemorySpace = dealii::MemorySpace::Host]’:
/home/ddong/Codes/PHiLiP/src/instantiate_vector_ad.cpp:26:1:   required 
from here
/home/ddong/Libraries/dealii/install/include/deal.II/lac/la_parallel_vector.templates.h:1651:7:
 
error: no matching function for call to 
‘std::complex::complex(const Sacado::Fad::DFad&)’
   AssertIsFinite(factor);
In file included from /usr/include/trilinos/Teuchos_ConfigDefs.hpp:94:0,
 from /usr/include/trilinos/Teuchos_PromotionTraits.hpp:45,
 from 
/usr/include/trilinos/Sacado_Fad_Exp_GeneralFadTraits.hpp:139,
 from /usr/include/trilinos/Sacado.hpp:52,
 from 
/home/ddong/Codes/PHiLiP/src/instantiate_vector_ad.cpp:1:
/usr/include/c++/7/complex:1512:3: note: candidate: constexpr 
std::complex::complex(const std::complex&)
   complex::complex(const complex& __z)
   ^~~
/usr/include/c++/7/complex:1512:3: note:   no known conversion for argument 
1 from ‘const Sacado::Fad::DFad’ to ‘const std::complex&’
/usr/include/c++/7/complex:1219:26: note: candidate: constexpr 
std::complex::complex(const std::complex&)
   _GLIBCXX_CONSTEXPR complex(const complex& __z)
  ^~~
/usr/include/c++/7/complex:1219:26: note:   no known conversion for 
argument 1 from ‘const Sacado::Fad::DFad’ to ‘const 
std::complex&’
/usr/include/c++/7/complex:1209:26: note: candidate: constexpr 
std::complex::complex(double, double)
   _GLIBCXX_CONSTEXPR complex(double __r = 0.0, double __i = 0.0)
  ^~~
/usr/include/c++/7/complex:1209:26: note:   no known conversion for 
argument 1 from ‘const Sacado::Fad::DFad’ to ‘double’
/usr/include/c++/7/complex:1207:26: note: candidate: constexpr 
std::complex::complex(std::complex::_ComplexT)
   _GLIBCXX_CONSTEXPR complex(_ComplexT __z) : _M_value(__z) { }
  ^~~
/usr/include/c++/7/complex:1207:26: note:   no known conversion for 
argument 1 from ‘const Sacado::Fad::DFad’ to 
‘std::complex::_ComplexT {aka __complex__ double}’
/usr/include/c++/7/complex:1202:12: note: candidate: constexpr 
std::complex::complex(const std::complex&)
 struct complex
^~~
/usr/include/c++/7/complex:1202:12: note:   no known conversion for 
argument 1 from ‘const Sacado::Fad::DFad’ to ‘const 
std::complex&’
/usr/include/c++/7/complex:1202:12: note: candidate: constexpr 
std::complex::complex(std::complex&&)
/usr/include/c++/7/complex:1202:12: note:   no known conversion for 
argument 1 from ‘const Sacado::Fad::DFad’ to 
‘std::complex&&’
In file included from 
/home/ddong/Libraries/dealii/install/include/deal.II/base/aligned_vector.h:22:0,
 from 
/home/ddong

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

[Praveen C]

If you are working only on Cartesian grids, these matrices will be same for all 
cells, except maybe for scalar factors. In this case, you can first compute 
them on a reference cell. For assembly, you loop over cells and perform a small 
matrix-vector product on each cell. Similar things can be done for face 
integrals.
best
praveen

> On 17-Sep-2019, at 12:33 PM, 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.
> 
> 
> 
> 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.
> 
> Before the any time steps, use a MeshWorker::loop each for each of the four 
> matrices: mass, differentiation, flux and boundary
> During each update
> Again use MeshWorker::loop, but this time only to calculate the normal 
> numerical flux.
> Use the normal numerical flux and the previous solution to obtain the RHS 
> using appropriate matrix-vector products
> Solve the system
> I have few question regarding this approach.
> Is it feasible
> Can it give significant improvement in performance over the case when 
> assembling is done for every time step
> (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/7fd51309-f137-421f-91bc-8c2d1c6b7ff1%40googlegroups.com
>  
> .
> 

-- 
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/8866E907-F297-4B65-AABD-D2F7E1BCE6D4%40gmail.com.

Re: [deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[Jean-Paul Pelteret]

H Doug and Wolfgang,

At the time when we introduced the AD framework I took a stab at quickly adding 
support for to our Vector class. To summarise, I hit the same issue and 
couldn’t find an elegant solution past to get past it (so hopefully this is the 
only outstanding problem to add support to Vector & friends). We have all of 
the tooling in place to convert AD numbers to something printable, but the 
issue is that we create circular inclusions in our headers when trying to use 
them, since the exceptions.h is a very fundamental header and number.h and 
almost all other headers require it.

My suggestion would be to not attack the problem in exceptions.h, but rather in 
vector.templates.h and la_parallel_vector.templates.h itself. So each call to 
AssertIsFinite(), for example over here 
,
 should have sent into it a number that’s pre-converted to a 
std::complex. You should be able to do this by invoking the conversion 
tool that will extract the value data for all AD numbers (so for free you’d 
include support for all Sacado and ADOL-C types):
AssertIsFinite(internal::NumberType>::value(s));

I’m pretty sure that will work. It’s a bit annoying, but I really think that 
might be the best way forward. If we go with this approach then I would suggest 
that we also add a note to the documentation of AssertIsFinite() to highlight 
this solution.

I hope this helps.

Best,
Jean-Paul

> On 18 Sep 2019, at 03:40, Doug  wrote:
> 
> so 'a' is of type Sacado::Fad::DFad. It then needs to call 
> 
>numbers::is_finite (Sacado::Fad::DFad), which doesn't exist. 
> It probably just tries to go through all of the overloads of 
> numbers::is_finite() and wants to see whether it can convert the 
> argument Sacado::Fad::DFad to the argument type of these 
> overloads. The error message you show then would just explain why this 
> one possibility (namely, numbers::is_finite(std::complex&)) 
> does not work. But that doesn't mean that this is the right overload 
> anyway -- I suspect that your compiler produces similar error messages 
> above or below the one you show for all of the other overloads, right? 
> 
> 
> True, the error message gets long pretty quickly, but the undefined is_finite 
> was another issue. Even if is_finite exists, the complex constructor is still 
> an issue.
>  
> I *think* that the solution is to simply provide an overload for 
>numbers::is_finite (const Sacado::Fad::DFad &x) 
> Can you try this? You could declare it in your own .cc file before you 
> #include  
> 
>  Attached is a tiny .cc file that simply instantiates the vector. Copy pasted 
> here here convenience.
> 
> #include 
> namespace dealii{
> namespace numbers{
> bool is_finite(const Sacado::Fad::DFad &x) {
> ¦   (void) x;
> ¦   return true;
> }   
> }}
> #include 
> #include 
> int main (int /*argc*/, char * /*argv*/[]){
> using namespace dealii;
> dealii::LinearAlgebra::distributed::Vector vector_double;
> using ADtype = Sacado::Fad::DFad;
> dealii::LinearAlgebra::distributed::Vector vector_ad;
> return 0;
> }
> 
> I also provide an function for numbers::is_finite (const 
> Sacado::Fad::DFad &x) to avoid the first set of errors. However, I 
> still get the error below
> 
> In file included from 
> /home/ddong/Libraries/dealii/install/include/deal.II/base/aligned_vector.h:22:0,
>  from 
> /home/ddong/Libraries/dealii/install/include/deal.II/lac/vector.h:22,
>  from 
> /home/ddong/Codes/PHiLiP/src/instantiate_vector_ad.cpp:13:
> /home/ddong/Libraries/dealii/install/include/deal.II/lac/la_parallel_vector.templates.h:
>  In instantiation of ‘dealii::LinearAlgebra::distributed::Vector MemorySpace>& dealii::LinearAlgebra::distributed::Vector MemorySpace>::operator*=(Number) [with Number = Sacado::Fad::DFad; 
> MemorySpace = dealii::MemorySpace::Host]’:
> /home/ddong/Codes/PHiLiP/src/instantiate_vector_ad.cpp:26:1:   required from 
> here
> /home/ddong/Libraries/dealii/install/include/deal.II/lac/la_parallel_vector.templates.h:1651:7:
>  error: no matching function for call to ‘std::complex::complex(const 
> Sacado::Fad::DFad&)’
>AssertIsFinite(factor);
> In file included from /usr/include/trilinos/Teuchos_ConfigDefs.hpp:94:0,
>  from /usr/include/trilinos/Teuchos_PromotionTraits.hpp:45,
>  from 
> /usr/include/trilinos/Sacado_Fad_Exp_GeneralFadTraits.hpp:139,
>  from /usr/include/trilinos/Sacado.hpp:52,
>  from 
> /home/ddong/Codes/PHiLiP/src/instantiate_vector_ad.cpp:1:
> /usr/include/c++/7/complex:1512:3: note: candidate: constexpr 
> std::complex::complex(const std::complex&)
>complex::complex(const complex& __z)
>^~~
> /usr/include/c++/7/complex:1512:3: note:   no known conversion for argument 1 
> from ‘const Sacado::Fad::DFad’ to ‘const std::complex&’
> /usr/include/c++/7/comp

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

[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.

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

[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)

[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)

[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.