Re: [deal.II] Re: Segmentation fault (core dumped) error

[Martin Kronbichler]

Hi Bruno,


2016-05-20 9:12 GMT-04:00 Ehsan :

#2  0x75cb614b in dealii::SparsityPattern::reinit
(this=this@entry=0x7fffb498, m=829992, n=829992, max_per_row=8032)
 at /deal_II/dealii-8.3.0/source/lac/sparsity_pattern.cc:251

It could be a memory problem. You try to create a matrix with 829,992
rows and up to 8,032 elements per row. This means 829,992*8,032 =
6,666,495,744 elements or 6,666,495,744 * 8 = 53,331,965,952 bytes =
53 Gb If you are only using one processor, you won't have enough
memory. However max_per_row seems excessively large, you can probably
find a sharper bound.


If I look at the current source code of SparsityPattern, I think there 
is indeed a bug: We define the variable max_vec_len to be of type 
size_type (= types::global_dof_index) whereas it actually should be 
std::size_t. We can have 32 bit integers but still more than 4G entries. 
I think that could be a problem here. Do you have time to have a look (I 
think one needs to check all the functions in SparsityPattern and the 
interface to SparseMatrix whether all types are correct)?


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: precomputation of mass matrix and stiffness matrix

[Martin Kronbichler]

Dear Marek,

Let's first look at the two distribute_local_to_global calls:


constraint_matrix_phase_constant.distribute_local_to_global(local_matrix,
dof_indices, system_matrix_phase_constant);



constraint_matrix_phase.distribute_local_to_global(local_matrix,
local_vector, dof_indices, jacobian_phase, 
residuum_phase);


I think that the problem is in the way you put together the entries for 
the constant part. When you assemble a single matrix, you most likely 
take the approach you use for 'jacobian_phase' with 
constraint_matrix_phase. Thus, in order to arrive at the same final 
system matrix, you need to use the same constraint matrix in both cases. 
The math behind this: In case of homogeneous constraints, the 
distribute_local_to_global method is a linear operation, so if the final 
local matrix is a sum of the two local matrices, you can equally well 
sum the two matrices.


So my assumption would be that the constraint_matrix_phase_constant 
object has different constraints than the other. This cannot work.


I assume you have different entries because you have some 
inhomogeneities on the latter which cannot be handled in the first call. 
I recommend to re-formulate the equations such that you only have 
homogeneous constraints, which then allows you to use the same 
constraint matrix on both parts of the matrix. This typically results in 
some additional right hand side terms including the left hand side terms.


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Step-35 Poiseuille flow in a 3D pipe

[Martin Kronbichler]

Dear Jiaqi,

Besides the literature suggested by Daniel, you should also read the paper

J.L. Guermond, P. Minev, J. Shen, An overview of projection methods for 
incompressible flows, Comput Methods Appl Mech Engrg 195, 2006.


Section 3.5 explains this question.

As a more general remark to your question I want to remind you that this 
forum is not primarily for asking questions about methods and their 
implications (go read the relevant papers) but rather more technical 
aspects regarding the implementation of finite element discretizations 
in deal.II. See also here: 
https://groups.google.com/forum/#!topic/dealii/GRZMUTLIm2I


Best,
Martin


On 03.08.2016 04:09, Jiaqi ZHANG wrote:

Hello David,

Thank you for the answer. I know Step-35 doesn't implement velocity 
correction, which


is still confusing me. Without correcting the velocity, it is not 
divergence free, does it make


any sense?

Best,
Jiaqi



2016-08-02 16:37 GMT-04:00 David Wells >:


Hi Jiaqi,

I have to admit that it has been about a year since I last really
looked at step-35 and I do not know the answer off the top of my
head for this. I saw that Daniel responded to this question
elsewhere; was that enough to figure things out?

Thanks,
David Wells


On Friday, July 29, 2016 at 9:19:29 PM UTC-4, Jiaqi ZHANG wrote:

Hello David Wells,

Sorry to bother you, but I have a problem about Step-35. I
have posted the question

a few days ago, and no one answered me, so I have to search in
the mailing list and found

that you have been using it a lot. I was wondering if you
could help me, and the following is

my question:

I read the references listed there and some other papers about
projection methods,

and know that after obtaining the intermediate values for
velocity, say v_n, we need to

correct the velocity by u_n = v_n - grad(2dt/3 phi_n), which
is mentioned in the introduction of Step-35,

but it seems that Step-35 only corrects the pressure, and I
cannot find anything in the code that implements the above
equation.

However, when I correct the velocity at every time step, the
algorithm turns out to be divergent (iterating more than 1000
times at some time step).

I am very confused, did I misunderstand the algorithm or the code?

Thanks in advance.


Best,
Jiaqi



在 2015年3月17日星期二 UTC-4下午8:14:21，David Wells写道：

Hi Joaquin,

I used step-35 a lot, so I may be able to assist you here.

For the velocity inflow, are you sure that what you have
will be zero flow at the edges? The XY slices don't look
quite circular to me, so this may be a problem if you
refine the mesh. Since you are using GMSH, it may be
easiest to use GMSH's extrusion abilities to extrude a
circle along the z axis to solve this.

I believe you handled the pressure correctly (as long as
the domain starts at z = 0 and ends at z = 10).

At least as of a few months ago, step-35 could not run in
3D due to a lack of a compressed sparsity pattern. I fixed
this in my own copy of step-35 (see
http://www.github.com/drwells/step-35a) (yeah, I know I
should submit a patch...). I tweaked the preconditioners
to improve performance, so my version might be useful to
look at.

Another thing you might want to look at is how Abner
implemented labels for the different boundaries. According
to my experiments (see

https://github.com/drwells/dealii-step-35a/blob/master/cylinderextruded.geo):
1 is the no-slip edges
2 is the inflow
3 is the outflow
4 is the cylinder (may not apply for you)

If your mesh file does not label things correctly then
things will get weird. I hope that this is helpful.

Thanks,
David Wells

On Tuesday, March 17, 2015 at 6:08:02 PM UTC-4, Joaquin
wrote:

Hi,

I am trying to run step-35 for the unsteady
incompressible viscous flow in a pipe, for 3D case.
The mesh file I used is attached. The graph of results
is attached. The expected results for any time should
be paraboloid (or parabolic in 2D), but, I can not get
that. The changes I did are:

*1. namespace EquationData:

*template 
double Velocity::value (const Point ,
 const unsigned int) const
{
  if (this->comp == 0)
   

Re: [deal.II] Re: Problem building dealii with PETSc

[Martin Kronbichler]

Dear Pete,

Let me just focus on two of your questions.

The first one is related to the OpenMPI installation on Ubuntu 16.04 
(and maybe other packages). You should not need more packages than what 
you list. However, it appears that the software shipped by Ubuntu is 
buggy. I can reproduce your error messages on my system 
(4.4.0-28-generic #47-Ubuntu SMP Fri Jun 24 10:09:13 UTC 2016 x86_64 
x86_64 x86_64 GNU/Linux):

$ cat test.cc
int main() {}
$ mpicxx -fuse-ld=gold test.cc
/usr/lib/openmpi/lib/libmpi_cxx.so: error: undefined reference to 
'opal_list_item_t_class'
/usr/lib/openmpi/lib/libmpi_cxx.so: error: undefined reference to 
'opal_class_initialize'

collect2: error: ld returned 1 exit status

This clearly is a bug in the OpenMPI implementation shipped by Ubuntu 
16.04. I think someone should open a bug on the Ubuntu (or Debian?) bug 
tracker but I lack the time to do so right now...


However, it only appears when using the gold linker. If one disables the 
gold linker, everything works fine.

$ mpicxx test.cc
$

Pete, can you try to go to the file 
"cmake/checks/check_01_compiler_features.cmake" in the deal.II source 
directory, search for the line

  ADD_FLAGS(DEAL_II_LINKER_FLAGS "-fuse-ld=gold")
(it should be at the very end of that file), comment it out, and try to 
build again with OpenMPI?


I think using the gold linker on this configuration is a bug in the MPI 
detection of deal.II (we really should use all compile and link flags). 
I have opened an issue here:

https://github.com/dealii/dealii/issues/2820

Regarding using MPICH on Ubuntu: I once tried to use that instead of 
OpenMPI but it is not really straight-forward. I think many other Debian 
packages rely on OpenMPI, so my experience is that you cannot really 
avoid to have both MPI packages installed. The only way for me to get 
MPICH to work rather than OpenMPI was to put symbolic links to MPICH in 
my PATH before the system paths, i.e., not really an elegant solution. 
(BTW: I cannot really recommend MPICH because my experience was that it 
got miserably slow if the MPI is over-allocated, i.e., more MPI 
processes than CPUs on my system. While you will not usually use this 
for production runs, I was not able to run the deal.II test suite in a 
reasonable time frame.)


Best,
Martin

On 13.07.2016 00:55, Pete Griffin wrote:

Wolfgang:

I installed:
sudo apt-get install libpetsc
sudo apt-get install petsc*
and
sudo apt-get install libopenlib
sudo apt-get install openlib*

I did not remove mpich since I assume since dealii was so specific in 
requiring the path info it shouldn't matter.


MPIhello world works ok.
gcc MPI-hello.c -I/usr/lib/openmpi/include 
/usr/lib/openmpi/lib/libmpi.so -o MPI-hello

./MPI-hello
Hello world from processor Ubuntu3, rank 0 out of 1 processors

I also have 7 questions relevant to future attempts at solving this 
problem. They are numbered below...


I ensured all the petsc and openmpi files as listed above were 
installed. I think they
probably were before to since the rest of the tests are the same.1) 
Are there any others

required that you know of?

Expand dealii-8.4.1.tar.gz to dealii-8.4.1-PETSc
cd /dealii-8.4.1-PETSc
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=/home/pgriffin/mylibs/dealii-8.4.1-PETSc 
-DDEAL_II_WITH_MPI=ON -DDEAL_II_WITH_PETSC=ON  ../

make -j8 install

This seemed to go O.K. as always for the non-compiled version. Screen 
config results are attached in DealiiConfigWithPETScOpenmpi.txt.


then

make test

I again as before got the following error (below) when running the 
dealii tests:


===< The dealii test error >

The results of the test are:

0% tests passed, 6 tests failed out of 6

Total Test time (real) =  92.47 sec

The following tests FAILED:
 1 - step.debug (Failed)
 2 - step.release (Failed)
 3 - affinity.debug (Failed)
 4 - mpi.debug (Failed)
 5 - tbb.debug (Failed)
 6 - step-petsc.debug (Failed)
Errors while running CTest

all of the errors are of the form:

CMakeFiles/Makefile2:3506: recipe for target 
'tests/quick_tests/CMakeFiles/step-petsc.debug.dir/all' failed
CMakeFiles/Makefile2:3545: recipe for target 
'tests/quick_tests/CMakeFiles/step-petsc.debug.run.dir/rule' failed


Complete screen results for test error in attached file 
DealiiTestWithPETScOpenmpi.txt.


2) Is this problem a dealii issue or should I contact petsc and/or 
openmpi development team, if mpich doesn't work.




I saw another person had the same error under UBUNTU 16.04 used mpich 
and was successful.


But as I said dealii refuses to find mpich and only sees openmpi.

I tried to override using:

cmake -DCMAKE_INSTALL_PREFIX=/home/pgriffin/mylibs/dealii-8.4.1-PETSc 
-DMPI_DIR=/usr/lib/mpich -DMPI_ARCH=x86_64 -DDEAL_II_WITH_MPI=ON 
-DDEAL_II_WITH_PETSC=ON  ../


3) Is this was the correct syntax to select an alternate MPI lib. My 
interpretation of the 

Re: [deal.II] Re: Problem building dealii with PETSc

[Martin Kronbichler]

Hi Pete,

I'm glad that it works. You may include my comments, the first line is 
the output of 'uname -srv' regarding my system information FYI.


Thanks!

Best,
Martin

On 13.07.2016 17:04, Pete Griffin wrote:

Thanks All!

Martin, after I commented out the line "ADD_FLAGS(DEAL_II_LINKER_FLAGS 
"-fuse-ld=gold")" in "cmake/checks/check_01_compiler_features.cmake", 
make test, worked. I ran step-17 and saw the speed improvement.

I will try (haven't done it before) to open a bug report on UBUNTU.

I will include your comments, as below, unless I hear from you otherwise.
I'm sure you have defined the problem better than I could!

===
... I can reproduce your error messages on my system (4.4.0-28-generic 
#47-Ubuntu SMP Fri Jun 24 10:09:13 UTC 2016 x86_64 x86_64 x86_64 
GNU/Linux):

$ cat test.cc
int main() {}
$ mpicxx -fuse-ld=gold test.cc
/usr/lib/openmpi/lib/libmpi_cxx.so: error: undefined reference to 
'opal_list_item_t_class'
/usr/lib/openmpi/lib/libmpi_cxx.so: error: undefined reference to 
'opal_class_initialize'

collect2: error: ld returned 1 exit status

This clearly is a bug in the OpenMPI implementation shipped by Ubuntu 
16.04. I think someone should open a bug on the Ubuntu (or Debian?) 
bug tracker but I lack the time to do so right now...


However, it only appears when using the gold linker. If one disables 
the gold linker, everything works fine.

$ mpicxx test.cc
$
=

Pete Griffin


On Sunday, July 10, 2016 at 10:54:54 PM UTC-4, Pete Griffin wrote:

I am trying to get PETSc running with dealii-8.4.1. (on Ubuntu
16.04) Reading the documentation it looks
like I don't need to have MPI. (???) I downloaded petsc-3.6.4 and
most recently compiled it with:

./configure --with-cc=gcc --with-cxx=g++ --with-fc=gfortran
--download-fblaslapack  --with-mpi=0

This resulted in the following being created in
petsc-3.6.4/arch-linux2-c-debug:

ls -R lib

lib:
libfblas.a  libflapack.a  libpetsc.so  libpetsc.so.3.6
 libpetsc.so.3.6.4  petsc  pkgconfig

lib/petsc:
conf

lib/petsc/conf:
configure.log  error.logfiles  make.log
 PETScConfig.cmake  petscvariables  reconfigure-arch-linux2-c-debug.py
configure.log.bkp  fblaslapack  gmake.log  modules petscrules
RDict.dbtest.log


lib/petsc/conf/modules:
petsc

lib/petsc/conf/modules/petsc:
3.6.4

lib/pkgconfig:
PETSc.pc

and

ls -R include/

include/:
mpi.modpetscdef.mod  petscdm.mod  
 petscksp.modpetscpcdef.modpetscsys.mod
petscaodef.mod petscdmcomposite.mod  petscfix.h  
petscmachineinfo.h  petscpc.mod   petsctsdef.mod
petscao.modpetscdmdadef.mod  petscisdef.mod  
petscmatdef.mod petscsnesdef.mod  petscts.mod
petscconf.hpetscdmda.mod petscis.mod  
 petscmat.modpetscsnes.mod petscvecdef.mod
petscconfiginfo.h  petscdmdef.mod  petsckspdef.mod  petsc.mod
  petscsysdef.mod petscvec.mod


Originally I compiled without the --with-mpi=0 and mpi libraries
were created. The MPI version
did not work! Since I am running on a single computer MPI was not
needed and tried without it.


=

The basic error is, with full details below:

-- Found PETSC
-- Could not find a sufficient PETSc installation: PETSc has to be
   configured with the same MPI configuration as deal.II.
-- DEAL_II_WITH_PETSC has unmet external dependencies.
CMake Error at cmake/configure/configure_3_petsc.cmake:110 (MESSAGE):

  Could not find the petsc library!

  Could not find a sufficient PETSc installation:

  PETSc has to be configured with the same MPI configuration as
deal.II, but
  found:

DEAL_II_WITH_MPI = OFF
PETSC_WITH_MPI   = (NOT FALSE)



I assume the -DDEAL_II_WITH_MPI=OFF switch would work since I
created PETsc with
  --with-mpi=0 and since mpi files were not created.



Thanks, Pete Griffin



=
THE RESULTS OF cmake for dealii were as follows:

=

$ cd dealii-8.4.1-PETSc/
/dealii-8.4.1-PETSc$ mkdir build
/dealii-8.4.1-PETSc$ cd build

=
/dealii-8.4.1-PETSc/build$ cmake
-DCMAKE_INSTALL_PREFIX=/home/pgriffin/mylibs/dealii-8.4.1-PETSc
-DPETSC_DIR=/petsc-3.6.4/arch-linux2-c-debug/lib
-DPETSC_ARCH=x86_64  

Re: [deal.II] Question concerning BlockSparsityPattern.copy_from() member function

[Martin Kronbichler]

Dear Dustin,

With regard to a lot of tutorial programms I would say that this is 
the usual way to work with /DynamicSparsityPatterns/, isn't it?
Yes, the description you give is exactly the way I would choose for the 
sparsity pattern.


So back to the question. Running the program the copy operation takes 
a huge amount of time. Just giving you some figures. With

ndofs_m = 823875,
ndofs_s = 1635075,
ndofs_bi = 25155,
meaning about 2.5 million dofs at all

it takes about 10 hours and 50 minutes to copy the 
/BlockDynamicSparsityPattern /in contrast to an assembling time of 2 
minutes and 38 seconds. So the question is if it is normal the copying 
operation taking so much time?
No, copying the sparsity pattern should in general be very fast. The 
functions that do this should be reasonably optimized so that you mostly 
pay for the memory access. It could be that we missed something for the 
block case, though. How is your computational setup, i.e., how many 
nonzero entries do you have in your matrix? Have you checked that you do 
not run out of memory and see a large swap time? How do the run times 
behave when you choose a smaller problem size? (I wonder if there is 
some higher than O(N) complexity somewhere.)



And also whether there is way to increase the performance?
By the way the program was compiled in release mode.
It should be possible to do the whole copy operation in say 2x the time 
it takes to zero a sparse matrix. If it's a bug we will fix it. It would 
be very helpful if you could provide us an example file that only 
contains the setup phase so we can investigate the issue further.


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Question concerning BlockSparsityPattern.copy_from() member function

[Martin Kronbichler]

Dear Dustin,



How is your computational setup, i.e., how many nonzero entries do
you have in your matrix?

I'm not sure if I understand what you mean. Do you mean the number of 
nonzero entries in /SparseMatrix/ or in the /BlockSparsityPattern/ or 
in the dynamic one? How can I get this information?
You can call DynamicSparsityPattern::n_nonzero_elements() to get the 
number of nonzero entries in the dynamic sparsity pattern. This method 
also exists in BlockSparsityPattern (or all sparsity patterns that 
inherit from BlockSparsityPatternBase):

https://dealii.org/developer/doxygen/deal.II/classBlockSparsityPatternBase.html

What I'm trying to understand here is what kind of properties your 
problem has - whether there are many nonzero entries per row and other 
special things that could explain your problems.


I just checked the 3D case of step-22 for the performance of 
BlockSparsityPattern::copy_from(BlockDynamicSparsityPattern) and the 
performance looks where I would expect it to be. It takes 1.19s to copy 
the sparsity pattern for a case with 1.6m DoFs (I have some 
modifications for the mesh compared to what you find online) on my 
laptop. Given that there are 275m nonzero entries in that matrix and I 
need to touch around 4.4 GB (= 4 x 275m x 4 bytes per unsigned int, once 
for clearing the data in the pattern, once for reading in the dynamic 
pattern, once for writing into the fixed pattern plus once for 
write-allocate on that last operation) of memory here, I reach 26% of 
the theoretical possible on this machine (~14 GB/s memory transfer per 
core). While I would know how to reach more than 80% of peak memory 
bandwidth here, this function is no way near being relevant in the 
global run time in any of my performance profiles. And I'm likely the 
deal.II person with most affinity to performance numbers.


Thus my interest in what is particular about your setup.


Have you checked that you do not run out of memory and see a large
swap time?

 I'm quiet sure that this is not the case/problem since I used one of 
our compute servers with 64 GB memory. Moreover, at the moment the 
program runs with an additional global refinement, i.e. about 16 
million dofs and only 33% of the memory is used. Swap isn't used at all.
That's good to know, so we can exclude the memory issue. Does your 
program use multithreading? It probably does in case you do not do 
anything special when configuring deal.II; the copy operation is not 
parallelized by threads but neither are almost all other initialization 
functions, so it should not become such a disproportionate timing here. 
10h for 2.5m dofs looks insane. I would expect something between 0.5 and 
10 seconds, depending on the number of nonzeros in those blocks.


Is there anything else special about your configuration or problem as 
compared to the cases presented in the tutorial? What deal.II version 
are you using, what is the finite element? Any special constraints on 
those systems?


Unfortunately this can not be done that easy. I have to reorganize 
things and kill a lot of superflous code. But besides that, I have a 
lot of other work to do. May be I can provide you an example file at 
the end of next week.
Let us know when you have a test case. I'm really curious what could 
cause this huge run time.


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Question concerning BlockSparsityPattern.copy_from() member function

[Martin Kronbichler]

Dear Dustin,

I had a look at your program on my machine with the most recent 
development version of deal.II. Here is the final timer output in 
release mode:


+-+++
| Total wallclock time elapsed since start| 1.02s ||
| |||
| Section | no. calls |  wall time | % of total |
+-+---+++
| Setup   | 1 | 0.618s |   60.% |
| copy_sparsity_pattern   | 1 | 0.0557s |   5.4% |
| setup_intial_fe_distribution| 1 | 0.0119s |   1.2% |
| setup_system| 1 | 0.606s |   59.% |
+-+---+++

As you can see, the copy operation is very cheap on my machine at 
0.0557s versus 20.2s on your machine. When looking at the profiles, I 
cannot see anything suspicious either. My primary suspecion right now is 
that you are on an old version of deal.II where we had bad code in that 
function. (But I don't remember having changed anything during the last 
year.) What version are you using?


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Segfault in EvaluatorTensorProduct::apply when AVX is enabled

[Martin Kronbichler]

Dear Stephen,

The problem is data alignment: You create an std::vector 
that internally arranges its data under the assumption that the start 
address of FEEvaluation is divisible by 32 (the length of the 
vectorization). If you put an FEEvaluation object on the stack, the 
compiler will automatically do it right. However, inside an std::vector 
it would be up to the std::vector to ensure this, but on usual x86-64 
machines it only aligns to 16 byte boundaries. This is also why SSE2 
works because it only needs 16-byte alignment.


The solution is to use AlignedVector scalar_vars instead of 
std::vector. The alternative is to wait until the pull request 3980 
(https://github.com/dealii/dealii/pull/3980) gets merged, grab the 
newest developer version because with that we will start using an 
external scratch data array that always has the correct alignment.


Best,
Martin


On 20.02.2017 20:50, Stephen DeWitt wrote:

Hello,
I recently realized that I should be using "-march=native" flag for 
optimal performance of matrix-free codes. The application that I've 
been using works fine with just SSE2, but with AVX enabled I'm getting 
a segfault. Step-48 works fine, so I don't think it is an installation 
issue.


The function where it occurs is similar to the "local_apply" function 
in step 48:

|
template
voidgetRHS(constMatrixFree,
   std::vector,
conststd::vector,
conststd::pair_range)const{

//initialize FEEvaulation objects
  std::vectorscalar_vars;

for(unsignedinti=0;i.

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Segfault in EvaluatorTensorProduct::apply when AVX is enabled

[Martin Kronbichler]

Dear Stephen,

The first problem you are seeing is a bug in AlignedVector::push_back.
We merged a pull request last night,
https://github.com/dealii/dealii/pull/3993, which should fix this issue.
The problem was that we incorrectly invoked placement new with default
constructor rather than the copy constructor.

Regarding your second approach: I am not sure what prevents this from
working. I have tried a similar approach and it does work in principle,
at least with the latest developer version. Apart from the fact that you
actually access invalid memory since you only call reserve(1) and never
extend the size to one.

Could you please try with the updated deal.II version and report back in
case you have further problems?

Thanks!

Best,
Martin


On 02/20/2017 11:01 PM, Stephen DeWitt wrote:
> Dear Martin,
> Thank you, that makes sense. 
>
> If I could trouble you with one more question -- I'm having a hard
> time using AlignedVector, and I can't find any examples of its use online.
>
> I'd like to build the AlignedVector with a series of push backs, but I
> get compiler errors if I try this:
>
> |
> dealii::AlignedVectorscalar_vars;
> typeScalar var1(data,0);
> scalar_vars.push_back(var1);
> |
>
> tellling me that:
> /Applications/deal.II.app/Contents/Resources/include/deal.II/base/aligned_vector.h:640:21:
> error: no matching constructor for initialization of
> 'dealii::FEEvaluation<2, 1, 2, 1, double>'
> new (_end_data) T;
>
> and:
> /Applications/deal.II.app/Contents/Resources/include/deal.II/base/aligned_vector.h:641:16:
> error: object of type 'dealii::FEEvaluation<2, 1, 2, 1, double>'
> cannot be assigned because its copy assignment operator is implicitly
> deleted
> *_end_data++ = in_data;
>
> I also tried:
>
> |
> dealii::AlignedVectorscalar_vars;
> scalar_vars1.reserve(1);
> typeScalar var1(data,0);
> scalar_vars1[0]=var1;
> |
>
> but still got the "error: object of type 'dealii::FEEvaluation<2, 1,
> 2, 1, double>' cannot be assigned because its copy assignment operator
> is implicitly deleted" error.
>
> How should I be making my AlignedVector?
>
> Thanks,
> Steve
>
>
>
> On Monday, February 20, 2017 at 3:23:09 PM UTC-5, Martin Kronbichler
> wrote:
>
> Dear Stephen,
>
> The problem is data alignment: You create an
> std::vector that internally arranges its data under
> the assumption that the start address of FEEvaluation is divisible
> by 32 (the length of the vectorization). If you put an
> FEEvaluation object on the stack, the compiler will automatically
> do it right. However, inside an std::vector it would be up to the
> std::vector to ensure this, but on usual x86-64 machines it only
> aligns to 16 byte boundaries. This is also why SSE2 works because
> it only needs 16-byte alignment.
>
> The solution is to use AlignedVector scalar_vars
> instead of std::vector. The alternative is to wait until the pull
> request 3980 (https://github.com/dealii/dealii/pull/3980
> <https://github.com/dealii/dealii/pull/3980>) gets merged, grab
> the newest developer version because with that we will start using
> an external scratch data array that always has the correct alignment.
>
> Best,
> Martin
>
>
> On 20.02.2017 20:50, Stephen DeWitt wrote:
>> Hello,
>> I recently realized that I should be using "-march=native" flag
>> for optimal performance of matrix-free codes. The application
>> that I've been using works fine with just SSE2, but with AVX
>> enabled I'm getting a segfault. Step-48 works fine, so I don't
>> think it is an installation issue.
>>
>> The function where it occurs is similar to the "local_apply"
>> function in step 48:
>> |
>> template
>> voidgetRHS(constMatrixFree<dim,double>,
>>  
>>  std::vector<dealii::parallel::distributed::Vector*>,
>>  
>>  conststd::vector<dealii::parallel::distributed::Vector*>,
>>conststd::pair<unsignedint,unsignedint>_range)const{
>>
>>   //initialize FEEvaulation objects
>>   std::vectorscalar_vars;
>>
>>   for(unsignedinti=0;i<num_var;i++){
>>   typeScalar var(data,i);
>>   scalar_vars.push_back(var);
>>   }
>>
>>   //loop over cells
>>  
>> for(unsignedintcell=cell_range.first;cell<cell_range.second;++cell){
>>
>>  // Initialize, read DOFs, and set evaulation flags  
>>  scalar_vars[varInfoListRHS[i].index].reinit(c

Re: [deal.II] FE_TraceQ using GaussLobatto nodes

[Martin Kronbichler]

Dear Praveen,

>
> When I use FE_TraceQ and am setting initial condition with
> VectorTools::interpolate, I get this error
>
> 
>
> An error occurred in line <127> of file
> 
> in function
>
> void dealii::VectorTools::interpolate(const Mapping
> &, const DoFHandlerType &, const Function typename VectorType::value_type> &, VectorType &, const
> dealii::ComponentMask &) [VectorType =
> dealii::TrilinosWrappers::MPI::Vector, dim = 2, spacedim = 2,
> DoFHandlerType = DoFHandler]
>
> The violated condition was: 
>
>
> (fe[fe_index].base_element(fe[fe_index].component_to_base_index(component_index).first).has_support_points())
> ||
> (fe[fe_index].base_element(fe[fe_index].component_to_base_index(component_index).first).dofs_per_cell
> == 0)
>
> The name and call sequence of the exception was:
>
> ExcNonInterpolatingFE()
>
> Additional Information: 
>
> You are attempting an operation that requires the finite element
> involved to be 'interpolating', i.e., it needs to have support points.
> The finite element you are using here does not appear to have those
> for the required components.
>
>
> which shouldn't occur since this is interpolating space.
You are right, we should allow for interpolation of initial condition.
We also allow this for FE_FaceQ. I will soon open a PR. On the other
hand, we will probably not be able to project onto FE_TraceQ with our
tools because this element does not define basis functions inside the
elements, only on the faces. (One can write manual project() functions
for FE_FaceQ, but I wouldn't know how to do this for FE_TraceQ.)

> Also in FE_TraceQ, there is a private member FE_Q fe_q.
> What is the use of this ?
It is strictly for internal use in order to simplify the implementation.
We use it for querying all kind of information, given that the degrees
of freedom of FE_TraceQ are a subset of the ones in FE_Q. The last
degrees of freedom of FE_Q, belonging to the interior of the element,
are left out, and polynomial spaces are interpreted as the restriction
of the FE_Q spaces to the surfaces.

Best,
Martin

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Question concerning BlockSparsityPattern.copy_from() member function

[Martin Kronbichler]

Hi Dustin,

Thanks for your file. This is indeed a bug in deal.II. Our algorithms 
for copying the sparsity pattern contain a part that scales 
quadratically in the number of rows in case most of the rows are empty, 
as is the case with some of your sparsity patterns. This is very bad if 
there are many rows. When using a patch, the run time in release mode of 
copy_sparsity_pattern goes from over 20s on my machine to 0.328s. A 
patch addressing this issue has been published on the deal.II site:


https://github.com/dealii/dealii/pull/3040

Best,
Martin


On 31.08.2016 13:20, Dustin Kumor wrote:

Dear all,

some additional information. Out of interest I switched back to the 
last release version of the deal.ii library 8.4.1 and made a run with 
this version. The problem is still there and even worse. I attached 
again the two log files for debug and release mode respectively.



Best
Dustin
--
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 
.

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] VectorTools::integrate_difference function

[Martin Kronbichler]

Oh, now I see: You should write "Solution", not "Solution<dim+1>".
The template parameters gives the space dimension of the function, not
the number of components. The latter is specified in the constructor of
your class Solution through the argument (dim+1) to Function.

Best,
Martin


On 09/20/2016 04:54 PM, JAEKWANG KIM wrote:
> thanks!! Full error message is as follow 
>
>
> */Users/jaekwangjk/repos/trial/sphere-deal.ii/mystokes.cc:1258:7:
> **error: **no*
>
> *  matching function for call to 'integrate_difference'*
>
>   VectorTools::integrate_difference (dof_handler,
>
> *  ^*
>
> */Users/jaekwangjk/repos/trial/sphere-deal.ii/mystokes.cc:1220:7:
> note: *in
>
>   instantiation of member function
>
>   'MyStokes::StokesProblem<2>::process_solution' requested here
>
>   process_solution ();
>
> *  ^*
>
> */Users/jaekwangjk/repos/trial/sphere-deal.ii/mystokes.cc:1298:21:
> note: *in
>
>   instantiation of member function 'MyStokes::StokesProblem<2>::run'
>
>   requested here
>
>flow_problem.run ();
>
> *^*
>
> */Users/jaekwangjk/Programs/dealii/include/deal.II/numerics/vector_tools.h:1940:8:
> note: *
>
>   candidate template ignored: could not match 'Function' against
> 'Solution'
>
>   void integrate_difference (const DoFHandler<dim,spacedim> ,
>
> *   ^*
>
> */Users/jaekwangjk/Programs/dealii/include/deal.II/numerics/vector_tools.h:1968:8:
> note: *
>
>   candidate template ignored: could not match 'dealii::hp::DoFHandler'
>
>   against 'dealii::DoFHandler'
>
>   void integrate_difference (const hp::DoFHandler<dim,spacedim> ,
>
> *   ^*
>
> */Users/jaekwangjk/Programs/dealii/include/deal.II/numerics/vector_tools.h:1925:8:
> note: *
>
>   candidate function template not viable: requires at least 7
> arguments, but
>
>   6 were provided
>
>   void integrate_difference (const Mapping<dim,spacedim>,
>
> *   ^*
>
> */Users/jaekwangjk/Programs/dealii/include/deal.II/numerics/vector_tools.h:1953:8:
> note: *
>
>   candidate function template not viable: requires at least 7
> arguments, but
>
>   6 were provided
>
>   void integrate_difference (const hp::MappingCollection<dim,spacedim>
> ,
>
>
>
>
> I will check the mapping as you advised for my other questions.!!
>
>
> 2016년 9월 20일 화요일 오전 9시 51분 13초 UTC-5, Martin Kronbichler 님의
> 말:
>
> Jaekwang,
>
> can you post the full error message? To me the call looks
> reasonable (this is a copy from step-7) and the function
> definition looks good as well.
>
> Ah, and in case you want to use this in the context of your other
> question regarding high order methods: Put a "mapping" as first
> argument to this call and make sure to use
> QGauss(fe.degree+2) or so, which increases the order of
> quadrature as the element order gets higher.
>
> Best,
> Martin
>
>
> On 09/20/2016 04:46 PM, JAEKWANG KIM wrote:
>>
>>  
>>
>> I wanted to evaluate Vectortool::integrate_difference to estimate
>> my error. 
>>
>>
>>
>> but I couldn't figure out what I need to enter for the third
>> component. (I marked it as red) 
>>
>>
>>
>>
>>  VectorTools::integrate_difference (dof_handler,
>>
>>  solution,
>>
>>  _Solution<dim+__1__>()_,
>>
>>  difference_per_cell,
>>
>>  QGauss(3),
>>
>>  VectorTools::L2_norm);
>>
>>
>> I have solution vector over my domain,(by the way I am solving
>> vector valued problems...)  
>>
>> I want to compare this with continuous exact function, "May be
>> Solution<dim+1>" , I have generated it as follow. 
>>
>> However, I got a error message that 
>>
>>
>> */Users/jaekwangjk/repos/trial/sphere-deal.ii/mystokes.cc:1258:7:
>> **error: **no*
>>
>> *  matching function for call to 'integrate_difference'*
>>
>> *
>> *
>>
>> What's might be wrong with my code..
>>
>>
>>  template
>>
>> classSolution : publicFunction
>>
>> {
>&

Re: [deal.II] CXX flags for SIMD and other optimization flags (GCC@6)

[Martin Kronbichler]

Hi Denis,

> From the matrix-free tutorials 37 and 48, I see that the recommended
> flags for VectorizedArrays with GCC are
>
> -DCMAKE_CXX_FLAGS="-march=native"
>
> How about using -O3, -ffast-math, -funroll-loops ? Any other
> recommended flags for GCC?
-march=native gives you AVX vectorization on most modern Intel CPUs
(starting from Sandy Bridge), which doubles the width of
VectorizedArray from 2 to 4. For computation bound algorithms
this almost doubles performance. Look for the line "-- Performing Test
DEAL_II_HAVE_AVX - Success" in the deal.II configuration to see whether
it gets enabled.

When comparing to a potential 2x speedup with AVX, -O3 helps only
little, and so does -ffast-math. In all my benchmark tests, their impact
has been on the level of noise. (One of my PhD students burns many
millions of CPU hours on big SandyBridge/Haswell clusters, which made me
spend hours on writing proposals, so we have checked.) -funroll-loops
should be enabled by default IIRC. For clang, I usually also enable
"-ffp-contract=fast" to enable fused multiply-add as it does not appear
per default. GCC does that in the default settings already it appears.

Best,
Martin

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Can we say that the higher order method, the more accurate?

[Martin Kronbichler]

Dear Jaekwank,

In addition to what Daniel and Wolfgang said: One does definitely benefit 
from going to higher degrees and deal.II is able to handle this (the 
accurate boundary representation by a mapping is one thing that is easily 
forgotten). In a recent preprint, we considered the flow around a cylinder, 
a standard benchmark test for the incompressible Navier-Stokes equations:
http://arxiv.org/pdf/1607.01323v1.pdf
Go to page 23, Table 1, to see how the solution quality increases as the 
polynomial degree is increased from 4 to 7.

Best,
Martin

On Sunday, September 18, 2016 at 5:57:16 PM UTC+2, JAEKWANG KIM wrote:
>
>
> Hello, I am a starter of dealii and am learning a lot these days with the 
> help of video lectures and tutorial examples. 
>
> I modified step-22 code (stokes flow code) into my own problem, the flow 
> around sphere.
>
> and I intend to evaluate the drag force (which is analytically given by 
> stokes equation) 
>
> My code reached quite close to the value since the absolute error  : 
> abs(drag_calculated-drag_exact)/drag_exact is around 10^(-3)
>
> However, I expected that if I input higher 'degree' I will receive more 
> accurate result, but it didn't
>
> Obviously Q2 is better than Q1. and Q3 is better than Q2. But Q4 or Q4 is 
> not better than Q2 or Q3? 
>
> Is there any reason on this? 
>
> (To be specific, if i say degree 2 , that mean I use (2+1) for velocity, 
> (2) for pressure, and (2+2) for Gauss integral
>
>
> Thank you 
>
> Jaekwang Kim  
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Can we say that the higher order method, the more accurate?

[Martin Kronbichler]

Dear Jaekwang,

if you are using the most recent developer version, you should
automatically get the Gauss-Lobatto version of the node distribution. We
made those point distributions the default for FE_Q(degree) this spring.

Regarding the limits of quadrature formulas: Depending on what exactly
you are doing, I would expect that the accuracy decreases as soon as you
go beyond degree 10 or so. I don't think that your problem is the
accuracy issues in terms of roundoff, because those issues would appear
first as you go beyond 1e-10.

Did you check the suggestions regarding the mapping? (I.e., you put a
"Mapping" argument to all the constructors of FEValues,
FEFaceValues, interpolate_boundary_values, etc.) Are the solver
tolerances tight enough? What happens if the mesh is refined? If this
still does not help, it might be good to share a small example where the
issue is observed.

Best,
Martin


On 09/20/2016 04:35 PM, JAEKWANG KIM wrote:
> thank for the reply!!
>
> my fe degree is declared as "fe (FE_Q(degree+1),
> dim,FE_Q(degree), 1)" and I used "QGauss  
> quadrature_formula(degree+2);" to calculate integral over the cell. 
>
> 2016년 9월 20일 화요일 오전 9시 26분 6초 UTC-5, Praveen C 님의 말:
>
>
> On Tue, Sep 20, 2016 at 7:49 PM, JAEKWANG KIM  > wrote:
>
> Can I ask more about "the limits of our implementation of our
> quadrature formulas?". 
> I wonder when it usually happens. 
>
>
> Once I calculated drag coefficient, with Q1, my error is 10%
> compare to exact solution. 
> However, I can 1.7% error when I use Q2 which is significant
> decrease!
> At Q3, I get 0.4% and At Q4 it starts to increase again 0.5%. 
> If I go higher, then my error is more than 100%. I really
> want to figure out why this happens
>
> To summarize 
> From Q1~Q3... it shows significant decrease in error
> but it is not anymore at Q4 and Q5 
>
>
> Hi
>
> If you are using FE_Q space with uniformly spaced support points,
> then there could be problem at higher degrees. Just to check, you
> should use Gauss nodes, e.g.
>
> FE_Q(QGaussLobatto<1>(degree+1))
>
> 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
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] VectorTools::integrate_difference function

[Martin Kronbichler]

Jaekwang,

can you post the full error message? To me the call looks reasonable
(this is a copy from step-7) and the function definition looks good as well.

Ah, and in case you want to use this in the context of your other
question regarding high order methods: Put a "mapping" as first argument
to this call and make sure to use QGauss(fe.degree+2) or so, which
increases the order of quadrature as the element order gets higher.

Best,
Martin


On 09/20/2016 04:46 PM, JAEKWANG KIM wrote:
>
>  
>
> I wanted to evaluate Vectortool::integrate_difference to estimate my
> error. 
>
>
>
> but I couldn't figure out what I need to enter for the third
> component. (I marked it as red) 
>
>
>
>
>  VectorTools::integrate_difference (dof_handler,
>
>  solution,
>
>  _Solution()_,
>
>  difference_per_cell,
>
>  QGauss(3),
>
>  VectorTools::L2_norm);
>
>
> I have solution vector over my domain,(by the way I am solving vector
> valued problems...)  
>
> I want to compare this with continuous exact function, "May be
> Solution" , I have generated it as follow. 
>
> However, I got a error message that 
>
>
> */Users/jaekwangjk/repos/trial/sphere-deal.ii/mystokes.cc:1258:7:
> **error: **no*
>
> *  matching function for call to 'integrate_difference'*
>
> *
> *
>
> What's might be wrong with my code..
>
>
>  template
>
> classSolution : publicFunction
>
> {
>
> public:
>
> 
>
> Solution () : Function(dim+1) {}
>
> 
>
> //const unsigned int b_id;
>
> virtualdoublevalue (constPoint   ,
>
>   constunsignedintcomponent = 0) const;
>
> 
>
> virtualvoidvector_value (constPoint ,
>
>Vector   ) const;
>
> 
>
> };
>
> 
>
> 
>
> template
>
> double
>
> Solution::value (constPoint  ,
>
> constunsignedintcomponent) const
>
> {
>
> Assert (component < this->n_components,
>
> ExcIndexRange (component, 0, this->n_components));
>
> 
>
> doubleU=-1.; //inflow velocity
>
> doubleobj_rad=0.2; //Object Radius
>
> doubler_sph = sqrt( pow(p[0],2) + pow(p[1],2) );
>
> doublephi = (-1) * acos(p[1]/r_sph);
>
> 
>
> doubleq_r_sph=U*cos(phi)*( 1+
> 0.5*pow(obj_rad,3)/pow(r_sph,3)-1.5*obj_rad/r_sph );
>
> doubleq_phi=(-1)*
> U*sin(phi)*(1-0.25*pow(obj_rad,3)/pow(r_sph,3)-0.75*obj_rad/r_sph);
>
> 
>
> doubleq_x=q_r_sph*sin(phi)+q_phi*cos(phi);
>
> doubleq_y=(-1)*q_r_sph*cos(phi)+q_phi*sin(phi);
>
> 
>
> 
>
> if(component == 0)
>
> returnq_x; 
>
> if(component == 1)
>
> returnq_y;
>
> else
>
> return0;  //
>
> 
>
> 
>
> }
>
>
> template
>
> void
>
> Solution::vector_value (constPoint ,
>
>Vector   ) const
>
> {
>
> for(unsignedintc=0; cn_components; ++c)
>
> values(c) = Solution::value (p, c);
>
> }
>
>
> -- 
> 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
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Gauss-Lobatto node distribution on FESYSTEM?

[Martin Kronbichler]

Dear Jaekwang,


> *FESystemfe; *
>
>
> and constructor mentions taylor hood space as follow 
>
>
>   StokesProblem::StokesProblem (constunsignedintdegree)
>
> :
>
> degree (degree),
>
>*fe (FE_Q(degree+**1), dim,*
>
> *FE_Q(degree), 1),*
>
You would write fe (FE_Q(QGaussLobatto<1>(degree+2)), dim,
FE_Q(QGaussLobatto<1>(degree+1)). In case you are using the
development version, you will automatically get the Gauss-Lobatto nodes
because we recently made them the default point distribution in higher
degree FE_Q elements. In case you are on 8.4 or older, the default
distribution is equidistant and you need the above constructor.

Best,
Martin

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Finite Difference with deal.ii and Parallel Distributed Triangulation

[Martin Kronbichler]

Dear Rajat,

You could do finite differences using the p4est wrapper in deal.II in
principle, but it would probably not be very efficient. You want to know
the indices which are not available in our wrapper at least, not sure
about p4est. To get them, I would build a DoFHandler based on some
element. FE_Q(1) would be the natural choice. Then you would need
to sort degrees of freedom in lexicographic way, which you can do
through a call to

DoFRenumbering::downstream(dof_handler, direction, true); 

where direction is something like (1, 1e-5, 1e-10). I'm not sure if we
ever tested this in parallel, but there is definitely a way for doing it.

More generally, using data structures like finite differences inside the
finite element capabilities of deal.II introduces quite some overhead
and you would lose most of the appealing features of finite differences.
Alternatively you might want to look into other, more FD-oriented
applications building on top of p4est, for example ForestClaw,
http://math.boisestate.edu/~calhoun/ForestClaw/

Best,
Martin


On 09/28/2016 06:30 AM, RAJAT ARORA wrote:
> Hello all,
>
> I am curious to know if it is feasible to use the P4est mesh structure
> with deal.ii to do a finite difference simulation.
>
> For a 2D structured grid, in order to do  finite differences, we need
> to know the index (or dof index) of vertices to the left as well  as
> to the right of a given grid point.
>
> Is there any function that can tell me the vertex index / dof index of
> the grid points which are just left and right to a chosen point. 
>
> Thanks.
> -- 
> 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
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] a few questions about matrix-free

[Martin Kronbichler]

Hi Denis,

I will put the description in the header of the file. PR underway. The
best documentation gets written when explaining things to others. :-)


> actually it’s not about components, it is about values at quadrature
> points for a vector of solution.
> Say I have two vectors representing scalar FE field and I want to have
> their values at quadrature points.
> That was my interpretation of the function below:
>
> void FEEvaluationBase
> <
> dim, n_components_, Number >::read_dof_values (   const 
> std::vector<
> VectorType > &src, 
>
>   
>   const unsigned int  first_index = |0| 
>
>   )
>
>
> since it takes a vector of VectorType.
> So the desired answer would be something like
> std::vector values = phi.get_value(q);
>
> Of course I could loop over the input vector and do
> phi.read_dof_values(src[i]) followed by phi.get_value(q) for each
> vector separately.
There are two possibilities for doing things. One is like you say, you
loop over the components in the input vector or create two FEEvaluation
objects and let them do similar things. Note though that you need to
call both 'phi.read_dof_values()' and 'phi.evaluate(true,false)' before
you can query the data on quadrature points.

The other possibility is to create a multi-component evaluator which was
along the lines of my last post. Say you have n_components vectors you
all want to read from, then you would set:
  FEEvaluation phi(matrix_free);

When you call 'phi.read_dof_values(src, 0);', it reads n_components
vectors one after another. The call to phi.get_value(q) then returns
Tensor<1,n_components,VectorizedArray > and you have all your
components sitting together.

Best,
Martin

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Alternative to local_range

[Martin Kronbichler]

Dear Praveen,

you can call IndexSet::nth_index_in_range(0) (in case 
IndexSet::n_elements() > 0) to get the first entry in the index set.


Best,
Martin


On 23.10.2016 07:38, Praveen C wrote:
I have an IndexSet of locally_owned_dofs and I want to know the 
starting index.


praveen

On Sun, Oct 23, 2016 at 9:56 AM, Wolfgang Bangerth 
> wrote:


On 10/22/2016 09:41 PM, Praveen C wrote:


local_range is marked deprecated


http://www.dealii.org/developer/doxygen/deal.II/classLinearAlgebra_1_1distributed_1_1Vector.html#ade4b13c22b96f11d6652bf2271602180



How can I get the same info without using this function ?


There should be a function that returns an IndexSet that contains
all of the elements locally stored.

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
.
For more options, visit https://groups.google.com/d/optout
.


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

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Interpolating p::d::Vector

[Martin Kronbichler]

Dear Praveen,

Called compress(VectorOperation::insert), but the element received 
from a remote processor, value 3.343515457745695e-15, does not match 
with the value -3.343515457745695e-15 on the owner processor 0



I dont understand why this error is triggered.
The error says that two processors have set the same vector entry to 
different values. On a VectorOperation::insert operation, the vector 
does not know which of the two conflicting values it should choose. This 
should not happen. Can you describe what kind of function you pass to 
the VectorTools::interpolate function? Given that the value is 
relatively small, I assume this happens because two processors decide to 
put the node location of a certain degree of freedom on a distance of 
machine epsilon apart. Given that they create different meshes, they 
will follow different paths down to a certain vertex along which they 
might refine differently.


The reason the check fails is that we measure a relative distance 
between the two entries 3.343515457745695e-15 and 
-3.343515457745695e-15, which is of course O(1) and fails. Since I'm 
hesitating to simply put an additional absolute check (what should 
happen if 1e-15 is a big number and small numbers are 1e-30?), I would 
try to get a global idea of the magnitude of expected entries in the 
vector by computing an infinity norm. This is a debug section only, so 
the cost is no problem. Any other idea?


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] LA::distributed::Vector initialized with MatrixFree and Kelly estimator

[Martin Kronbichler]

Hi Denis,

I wonder what is the best workaround given the fact that i have quite 
a number of such vectors?
You only see this for KellyErrorEstimator, right? Given that the error 
estimation is much more expensive than the operator evaluation anyway, I 
would suggest you do the copy exactly the way you describe here:


Probably the only way is to copy content into a temporary vector, then 
call reinit(locally_owned, locally_relevant, mpi) and assign the 
content back.

Then i would use such vectors in Kelly and solution transfer.

You would of course also invoke vector.update_ghost_values().

Any better ways to accomplish this? Maybe implement/add 
LA::distributed::Vector::set_ghost_set(const IndexSet &) but not sure 
it's worth the effort.
I think the current interface is pretty good after all. Doing one copy 
is not that bad. Trilinos does copy the vector to be ghosted for the 
parallel sparse matrix-vector product in every operator application. If 
you want to change the ghost set, you need to throw away the vector as 
you note in the next post. Otherwise, one cannot change the ghost range 
of the vector. In addition, changing just the ghost side inside the 
vector is extremely dangerous. It took me a few iterations to get the 
`partitioner_is_compatible` check work somewhat reliably. The main 
problem is that local index spaces with more than one set of ghosts are 
ambiguous per definition. Besides deal.II I know two more big projects 
that have local and global MPI indices and everyone does the same 
struggles: You want local indices for performance, but it's so easy to 
get things wrong.


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: LA::distributed::Vector initialized with MatrixFree and Kelly estimator

[Martin Kronbichler]

Wouldn't it be easier to just create a copy of the vector by

ghosted_vector.reinit(locally_owned_set, locally_relevant_set, comm);
ghosted_vector = solution;
ghosted_vector.update_ghost_values();

not in my case, because i have  a lot of vectors and if I can avoid, i would 
prefer not to create a ghost equivalent for each one of them.
But you are right, if there is only a single solution vector, your approach 
makes more sense and is also consistent with what we normally do
with Trilinos/Petsc, i.e. step-40.
That makes sense, as does your approach with changing the ghost set of 
the vectors in case you don't want to spent two solution vectors. As 
long as you don't reuse those vectors in a matrix free loop again you're 
fine.


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: LA::distributed::Vector initialized with MatrixFree and Kelly estimator

[Martin Kronbichler]

Hi Denis,

Wouldn't it be easier to just create a copy of the vector by

ghosted_vector.reinit(locally_owned_set, locally_relevant_set, comm);
ghosted_vector = solution;
ghosted_vector.update_ghost_values();

The reason why there is a VectorView in the MatrixFreeOperators::Base 
class is that I want to modify the original vector which I don't assume 
you need in your code. (It would actually be dangerous if you used a 
vector with the locally relevant dofs as ghosts because that gets the 
local indices wrong as compared to the storage of MatrixFree-compatible 
vectors.)


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] trouble running new step-55

[Martin Kronbichler]

Dear Franck,

I justed checked on my machine and it works both on an installation that 
only provides Trilinos and on an installation that has both (and where 
PETSc is selected). It looks like your installation is 8.4.0 whereas I 
think that step-55 requires the current developer version 8.5.0-pre. The 
compile error complains about missing classes in the LA::MPI namespaces 
which might be related to some changes we have done recently. Can you 
check with the developer version?


Best,
Martin


On 17.10.2016 12:12, 'franck75' via deal.II User Group wrote:

I am having trouble when trying to  run new step-55.

I am having a very long error message. Here is one piece (see attached 
file).


Does anyone able to run it?

Best
Franck


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

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Find dofs supported on each face and quadrature point

[Martin Kronbichler]

Dear Praveen,

No, we do not have information about these properties available in the 
library and it cannot be easily introduced either. The reason is that 
the particular "spectral" evaluation of polynomials with roots in the 
quadrature points depends on the interaction of the finite element and 
the quadrature formula, which are separated in deal.II and only come 
together for FEValues.


The way you describe in is almost exactly the way I would do it. At 
first sight, there only appears to be a minor bug in your code:




TableIndices<3> size_u (GeometryInfo::faces_per_cell,

  fe_u.degree+1,

  fe_u.degree+1);

face_dof_index_u.reinit (size_u);

In the code below, you use the second index for the quadrature index on 
the face, whose size is (fe_u.degree+1)^(dim-1) or simply 
quadrature_formula.size() if you move up that part a few lines.


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Workstream, number threads

[Martin Kronbichler]

Hi Praveen,

Do you happen to use a BLAS that uses threading? I recommend switching 
to non-threaded BLAS (or "export OMP_NUM_THREADS=1"). deal.II applies 
threading on a higher level and a threaded BLAS is counter-productive in 
98% of cases.


Best,
Martin


On 29.12.2016 11:30, Praveen C wrote:

Hello Daniel

I recompiled dealii without threads.

When I run step-40 single process

./step-40

I see upto 2000% cpu usage during cycle 7 when viewed in top.

Best
praveen

On Thu, Dec 29, 2016 at 3:51 PM, Daniel Arndt 
> wrote:


Praveen,

unsignedintn_threads = 1;

Utilities::MPI::MPI_InitFinalize mpi_initialization (argc,
argv, n_threads);


But when I run this code, each process is using 1700% of CPU.
The code does not progress in its iterations also. What am I
doing wrong ?

This looks correct. Can you verify that you don't have that
problem when you are running step-9 (for WorkStream) and step-40
(for MPI)?

Do you have a minimal working example that shows this problem?

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
.
For more options, visit https://groups.google.com/d/optout
.


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

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Internal instability of the GMRES Solver / Trilinos

[Martin Kronbichler]

Dear Pascal,

No, you do not need to try the other solution. I'm glad I could help.
(This asserts the approach that we need to be careful with the vector
pool between different calls.)

Best,
Martin


On 16.03.2017 15:21, Pascal Kraft wrote:
> Hi Martin,
>
>  I have tried a version
> with 
> GrowingVectorMemory::release_unused_memory()
> at the end of each step and removed my change to trilinos_vector.cc
> l.247 (back to the version from dealii source) and it seems to work
> fine. I have not tried the other solution you proposed, should I?
> Would the result help you?
>
> Thank you a lot for your support! This had been driving me crazy :)
>
> Best,
> Pascal
>
> Am Donnerstag, 16. März 2017 08:58:53 UTC+1 schrieb Martin Kronbichler:
>
> Dear Pascal,
>
> You are right, in your case one needs to call
> 
> GrowingVectorMemory::release_unused_memory()
> rather than for the vector. Can you try that as well?
>
> The problem appears to be that the call to SameAs returns
> different results for different processors, which it should not
> be, which is why I suspect that there might be some stale
> communicator object around. Another indication for that assumption
> is that you get stuck in the initialization of the temporary
> vectors of the GMRES solver, which is exactly this kind of situation.
>
> As to the particular patch I referred to: It does release some
> memory that might have stale information but it also changes some
> of the call structures slightly. Could you try to change the
> following:
>
> if(vector->Map().SameAs(v.vector->Map()) == false)
>
> to
>
> if(v.vector->Map().SameAs(vector->Map()) == false)
>
> Best, Martin
>
> On 16.03.2017 01:28, Pascal Kraft wrote:
>> Hi Martin,
>> that didn't solve my problem. What I have done in the meantime is
>> replace the check in line 247 of trilinos_vector.cc with true. I
>> don't know if this causes memory leaks or anything but my code
>> seems to be working fine with that change. 
>> To your suggestion: Would I have also had to call the templated
>> version for BlockVectors or only for Vectors? I only tried the
>> latter. Would I have had to also apply some patch to my dealii
>> library for it to work or is the patch you talked about simply
>> that you included the functionality of the call
>> 
>> GrowingVectorMemory::release_unused_memory()
>> in some places?
>> I have also wanted to try MPICH instead of OpenMPI because of a
>> post about an internal error in OpenMPI and one of the functions
>> appearing in the call stacks sometimes not blocking properly.
>>     Thank you for your time and your fast responses - the whole
>> library and the people developing it and making it available are
>> simply awesome ;)
>> Pascal
>> Am Mittwoch, 15. März 2017 17:26:23 UTC+1 schrieb Martin Kronbichler:
>>
>> Dear Pascal,
>>
>> This problem seems related to a problem we recently worked
>> around in https://github.com/dealii/dealii/pull/4043
>> <https://github.com/dealii/dealii/pull/4043>
>>
>> Can you check what happens if you call
>> 
>> GrowingVectorMemory::release_unused_memory()
>>
>> between your optimization steps? If a communicator gets stack
>> in those places it is likely a stale object somewhere that we
>> fail to work around for some reason.
>>
>> Best, Martin
>>
>> On 15.03.2017 14:10, Pascal Kraft wrote:
>>> Dear Timo,
>>> I have done some more digging and found out the following.
>>> The problems seem to happen in trilinos_vector.cc between
>>> the lines 240 and 270.
>>> What I see on the call stacks is, that one process reaches
>>> line 261 ( ierr = vector->GlobalAssemble (last_action); )
>>> and then waits inside this call at an MPI_Barrier with the
>>> following stack:
>>> 20  7fffd4d18f56
>>> 19 opal_progress()  7fffdc56dfca
>>> 18 ompi_request_default_wait_all()  7fffddd54b15
>>> 17 ompi_coll_tuned_barrier_intra_recursivedoubling()
>>>  7fffcf9abb5d
>>> 16 PMPI_Barrier()  7fffddd68a9c
>>> 15 Epetra_MpiDistributor::DoPosts()  7fffe4088b4f
>>> 14 Epetra_MpiDistributor::Do()  7fffe4089773
>>> 13 Epetra_DistObject::DoTransfer()  7fffe400a96a
&g

Re: [deal.II] Re: Internal instability of the GMRES Solver / Trilinos

[Martin Kronbichler]

Dear Pascal,

This problem seems related to a problem we recently worked around in 
https://github.com/dealii/dealii/pull/4043


Can you check what happens if you call 
GrowingVectorMemory::release_unused_memory()


between your optimization steps? If a communicator gets stack in those 
places it is likely a stale object somewhere that we fail to work around 
for some reason.


Best,
Martin


On 15.03.2017 14:10, Pascal Kraft wrote:

Dear Timo,

I have done some more digging and found out the following. The 
problems seem to happen in trilinos_vector.cc between the lines 240 
and 270.
What I see on the call stacks is, that one process reaches line 261 
( ierr = vector->GlobalAssemble (last_action); ) and then waits inside 
this call at an MPI_Barrier with the following stack:

20  7fffd4d18f56
19 opal_progress()  7fffdc56dfca
18 ompi_request_default_wait_all()  7fffddd54b15
17 ompi_coll_tuned_barrier_intra_recursivedoubling()  7fffcf9abb5d
16 PMPI_Barrier()  7fffddd68a9c
15 Epetra_MpiDistributor::DoPosts()  7fffe4088b4f
14 Epetra_MpiDistributor::Do()  7fffe4089773
13 Epetra_DistObject::DoTransfer()  7fffe400a96a
12 Epetra_DistObject::Export()  7fffe400b7b7
11 int Epetra_FEVector::GlobalAssemble()  7fffe4023d7f
10 Epetra_FEVector::GlobalAssemble()  7fffe40228e3

The other (in my case three) processes are stuck in the head of the 
if/else-f statement leading up to this point, namely in the line
if(vector->Map().SameAs(v.vector 
->Map()) 
== false)

inside the call to SameAs(...) with stacks like
15 opal_progress() 7fffdc56dfbc 14 ompi_request_default_wait_all() 
7fffddd54b15 13 ompi_coll_tuned_allreduce_intra_recursivedoubling() 
7fffcf9a4913 12 PMPI_Allreduce() 7fffddd6587f 11 
Epetra_MpiComm::MinAll() 7fffe408739e 10 Epetra_BlockMap::SameAs() 
7fffe3fb9d74
Maybe this helps. Producing a smaller example will likely not be 
possible in the coming two weeks but if there are no solutions until 
then I can try.

Greetings,
Pascal
--
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 
.

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: CMake extra libraries?

[Martin Kronbichler]

This code triggers a safety mechanism that I introduced recently (April 
26). We also see this on one of the test machines: 
https://cdash.kyomu.43-1.org/viewBuildError.php?buildid=8758
but I previously thought this is only due to a leftover in the 
configuration and I did not see it in real life before.


Victor, are you on a clean deal.II build directory without previous content?

Could you also post the last 20 lines or so regarding the instruction 
extension in the output of the intel compiler you are using with


|icpc |||-xCORE-AVX2 -axMIC-AVX512,CORE-AVX512| -dM -E -x c++ /dev/null |

(see also the discussion we had a short while ago here: 
https://github.com/dealii/dealii/pull/4440 )


Best,
Martin

On 07.07.2017 14:21, Bruno Turcksin wrote:

This is very strange. Do you set the AVX flags yourself?

@Martin have you ever seen something like that?

Bruno

2017-07-07 3:50 GMT-04:00 Victor Eijkhout :

Hm. That second compile line should be:









[ 63%] Building CXX object
source/numerics/CMakeFiles/obj_numerics_debug.dir/vector_tools_boundary.cc.o


cd /tmp/dealii-build/source/numerics &&
/opt/intel/compilers_and_libraries_2017.4.196/linux/mpi/intel64/bin/mpicxx
-DDEBUG -DTBB_DO_ASSERT=1 -DTBB_USE_DEBUG
-I/tmp/dealii-build/source/numerics -I/tmp/dealii-build/include -I/admi\


n/build/rpms/BUILD/dealii-git20170615/include
-I/admin/build/rpms/BUILD/dealii-git20170615/bundled/tbb41_20130401oss/include
-I/admin/build/rpms/BUILD/dealii-git20170615/bundled/umfpack/UMFPACK/Include
-I/admin/build/rpms/BUILD/d\


ealii-git20170615/bundled/umfpack/AMD/Include
-I/admin/build/rpms/BUILD/dealii-git20170615/bundled/muparser_v2_2_4/include
-I/opt/intel/compilers_and_libraries_2017.4.196/linux/mpi/intel64/include
-I/home1/apps/intel17/impi17_0/p\


etsc/3.7/knightslanding/include
-I/home1/apps/intel17/impi17_0/trilinos/12.10.1/include
-I/opt/apps/intel17/boost/1.64/include
-I/opt/apps/intel17/impi17_0/parallel-netcdf/4.3.3.1/x86_64/include
-I/opt/apps/intel17/impi17_0/phdf5\


/1.8.16/x86_64/include
-I/opt/intel/compilers_and_libraries/linux/mkl/include
-I/home1/apps/intel17/impi17_0/petsc/3.7/include
-I/home1/apps/intel17/impi17_0/p4est/2.0/include
-I/home1/apps/intel17/impi17_0/slepc/3.7/knightslandi\


ng/include -I/home1/apps/intel17/impi17_0/slepc/3.7/include  -fpic -ansi -w2
-wd21 -wd68 -wd135 -wd175 -wd177 -wd191 -wd193 -wd279 -wd327 -wd383 -wd981
-wd1418 -wd1478 -wd1572 -wd2259 -wd2536 -wd3415 -wd15531 -wd111 -wd128
-wd185\


  -wd280 -qopenmp-simd -std=c++14 -std=c++14 -Wno-return-type
-Wno-parentheses -O0 -g -gdwarf-2 -grecord-gcc-switches -std=c++14 -g -O0 -o
CMakeFiles/obj_numerics_debug.dir/vector_tools_boundary.cc.o -c
/admin/build/rpms/BUILD/dea\


lii-git20170615/source/numerics/vector_tools_boundary.cc


In file included from
/admin/build/rpms/BUILD/dealii-git20170615/include/deal.II/matrix_free/matrix_free.h(24),


  from
/admin/build/rpms/BUILD/dealii-git20170615/include/deal.II/numerics/vector_tools.templates.h(43),


  from
/admin/build/rpms/BUILD/dealii-git20170615/source/numerics/vector_tools_boundary.cc(17):


/admin/build/rpms/BUILD/dealii-git20170615/include/deal.II/base/vectorization.h(46):
error: #error directive: "Mismatch in vectorization capabilities: AVX was
detected during configuration of deal.II and switched on, but it is ap\


parently not available for the file you are trying to compile at the moment.
Check compilation flags controlling the instruction set, such as
-march=native."


   #error "Mismatch in vectorization capabilities: AVX was detected during
configuration of deal.II and switched on, but it is apparently not available
for the file you are trying to compile at the moment. Check compilation
flags \


controlling the instruction set, such as -march=native."



--
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 a topic in the
Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit
https://groups.google.com/d/topic/dealii/msEQFzuoxf0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to
dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] segmentation fault when calling make_vectorized_array

[Martin Kronbichler]

Hi Denis,




no, it's just a member variable VectorizedArray:

VectorizedArray a;


But then the class/struct that holds this member variable is not 
properly aligned, you probably have it inside an std::vector or 
allocated it with 'new'. You will need to put it into an AlignedVector 
or similar type. Alternatively you need to wrap the member that is of 
type VectorizedArray into an AlignedVector.




|
  >│0x773913cb ::Operator()+35>   vmovaps %ymm0,0xe0(%rbx)

|


This is the offending instruction. The compiler tries to store to the 
address given by an offset of 0xe0 (224 bytes) to the value in %rbx, 
using an 'aligned' store operation (vmovAps) that assumes 32 byte 
alignment. If you print the value of %rbx (in gdb with 'print $rbx'), 
I'm pretty sure you get an address that is divisible by 16 but not by 
32. It is a pity that the compiler / API of the intrinsics forces the 
compiler to use an aligned access instruction whenever there is a 
variable of type VectorizedArray but there's nothing we can change and 
we have to live with that...


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] segmentation fault when calling make_vectorized_array

[Martin Kronbichler]

Hi Denis,

what is a? Is it an element of an array? If yes, is the array of type 
AlignedVector (or derived from that).


You get this kind of error when there is a load or store to an address 
that is not aligned by the size of the vectorized array, 32 bytes in 
your case. The data types VectorizedArray make the compiler insert 
'aligned load/store' type of operations, which cannot always be guaranteed.


Can you try to look at the problem in a debugger, e.g. gdb and post the 
offending assembler line and a few before that to see the context?


Best,
Martin


On 12.07.2017 17:41, Denis Davydov wrote:

Hi all,

With the current master (182e974) I get a weird segmentation fault for 
*float* numbers when calling


|
 a =make_vectorized_array(0.);
|

inside a constructor of my class where

|
VectorizedArraya;
|

I see this only on a cluster (Intel Xeon Ivy Bridge), deal.II is 
compiled with *-march=native* with GCC 4.8.5 in Debug mode (obviously 
segmentation fault is there in Release as well).


Unfortunately, I wiped the old (working) installation and can't check 
it now,

but I think this used to work with *-march=native* a few months back.

Relevant config tests are:

|
--PerformingTestDEAL_II_HAVE_SSE2
--PerformingTestDEAL_II_HAVE_SSE2 -Success
--PerformingTestDEAL_II_HAVE_AVX
--PerformingTestDEAL_II_HAVE_AVX -Success
--PerformingTestDEAL_II_HAVE_AVX512
--PerformingTestDEAL_II_HAVE_AVX512 -Failed
--PerformingTestDEAL_II_HAVE_OPENMP_SIMD
--PerformingTestDEAL_II_HAVE_OPENMP_SIMD -Failed
|

If I try to replicate what *make_vectorized_array* does, then it works 
fine whereas calling it gives segfault:


|
constvalue_type u =0.;
VectorizedArrayresult;
  result =u;
// ^^ copy paste of make_vectorized_array
  a =u;
VectorizedArraydummy =make_vectorized_array(0.);// <- 
quick test for doubles
  a =make_vectorized_array(0.);<=Segmentationfault 
here!

|

I am puzzled with this... Any ideas where to dig?

Cheers,
Denis
--
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 
.

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] segmentation fault when calling make_vectorized_array

[Martin Kronbichler]

Hi Denis,


> While on the topic of alignment, do I need to be more careful about using
>
> template
> structQPData
> {
>Table<2, VectorizedArray > values;
> }
> QPData qp_data;
>
> as a member variable in my class?

Table is fine because it is based on AlignedVector,
so it will use the correct alignment internally.

Best,
Martin

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Refining a hyper_ball with a spherical_manifold produces nan values

[Martin Kronbichler]

Dear Alex,

you cannot set a SphericalManifold for the innermost part of a 
hyper_ball because the transformation degenerates for radius=0. You need 
to restrict the manifold to only the boundary or only the outer part of 
the hyper_ball.


Best,
Martin


P.S. for future reference, if I have an issue formulated like this, 
should I raise it on the deal.II repo or here in the mailing list?

The mailing list is the correct place for this kind of questions.

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Can we say that the higher order method, the more accurate?

[Martin Kronbichler]

Dear Howe,

How did you run your simulation? From your picture, it appears that a 
higher order method is worse at higher degrees than a lower order 
method, which does not match with my experience. If that were the case, 
nobody would use high orders. However, you need to bring many pieces in 
place to really get to the benefit of the high order method for somewhat 
more complicated examples such as the flow around a cylinder. Here is a 
list of things to look at:


- Do you use a high-order polynomial mapping MappingQ of the same or 
higher degree as the interpolation space? Do you use this mapping in all 
routines that evaluate quantities, such as the usual assembly, the 
computation of the lift/drag, and so on?
- Do you use a manifold description that extends into the domain? (Look 
into TransfiniteInterpolationManifold.) Without, you will not get more 
than third order convergence.
- Do you have a good mesh around the area of interest? Flows around 
cylinders tend to be really really sensitive to the mesh quality around 
the cylinder.


For the Navier-Stokes equations around the cylinder, if everything is 
done right one gets significantly improved results in terms of accuracy 
over the number of degrees of freedom up to degree (6,5) 
(velocity,pressure). Beyond that picture is less clear. At least with 
the meshes that we tried in our group it was not worth to go beyond. You 
can have a look a our results in section 5.4 and Figs. 9 and 10 of this 
preprint:

https://arxiv.org/pdf/1706.09252.pdf

Best,
Martin


On 09.08.2017 09:01, Howe wrote:

Dear Jaekwang

Have you solved this problem? If yes, Could pls share your solution 
with us?
I am simulating a steady state flow over a cylinder, and the drag/lift 
coefficient shows an unexpected trend of change as i increase the 
discretization order and refine the mesh.




As is shown in the figure, the Cd increased as the cells increased for 
all the discretization orders, however, for a fixed cells, the Cd 
decreased as the discretization order increased.


In my opinion, to increase the order and refine the mesh should both 
make the approximation more close to the exact solution, thus should 
have the same trend of change.



在 2016年9月18日星期日 UTC+8下午11:57:16，Jaekwang Kim写道：


Hello, I am a starter of dealii and am learning a lot these days
with the help of video lectures and tutorial examples.

I modified step-22 code (stokes flow code) into my own problem,
the flow around sphere.

and I intend to evaluate the drag force (which is analytically
given by stokes equation)

My code reached quite close to the value since the absolute error
 : abs(drag_calculated-drag_exact)/drag_exact is around 10^(-3)

However, I expected that if I input higher 'degree' I will receive
more accurate result, but it didn't

Obviously Q2 is better than Q1. and Q3 is better than Q2. But Q4
or Q4 is not better than Q2 or Q3?

Is there any reason on this?

(To be specific, if i say degree 2 , that mean I use (2+1) for
velocity, (2) for pressure, and (2+2) for Gauss integral


Thank you

Jaekwang Kim

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

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Can we say that the higher order method, the more accurate?

[Martin Kronbichler]

Dear Howe,

Regarding 1., the main question is whether you added the 'mapping' 
variable as the first arguments to all FEValues and FEFaceValues 
constructors, to calls to VectorTools::interpolate* and 
VectorTools::integrate_difference calls and the like, i.e., all routines 
that internally construct an FEValues or FEFaceValues object. There is 
defaults for many of those data structures using linear mappings, 
MappingQ1, but you don't want to use them but rather the high order 
description.


Regarding 2, I guess you are using a curved boundary description as 
explained in step-1? You need to tell the triangulation to use a curved 
description.


Regarding 3, your mesh looks quite good. You would probably want to use 
a volume manifold on the whole circular part of the domain, though. And 
ideally a TransfiniteInterpolationManifold on the regions where you go 
from the circle to the straight ends, but the latter is not yet 
available in any of the releases yet, only the github code, and it is 
not critical and you should get better results up to degree 3 at least.


Best,
Martin


On 09.08.2017 10:06, Howe wrote:

Dear Martin,

Thanks for your rapid response.
1.
The MappingQ is set to be the same as the order of velocity, as is 
shown in the following code snippet:


/template 
//NS::NS (ParameterHandler )
//:
//   parameters (),
//   degree (prm.get_integer("pressure degree")),
//   fe( FE_Q(QGaussLobatto<1>(degree+2)), dim,
// FE_Q(QGaussLobatto<1>(degree+1)), 1),
//   fe_scalar (FE_Q(QGaussLobatto<1>(degree+2))),
//   dof_handler (triangulation),
//   dof_handler_scalar (triangulation),
//mapping (degree+2),
//   computing_timer (std::cout,
//   TimerOutput::summary,
//   TimerOutput::wall_times)/


I am not quite sure whether the computation of the lift/drag in my 
code is right, and my implementation  is almost the same as the one in 
this post: 
https://groups.google.com/d/msg/dealii/rS6soTb69ig/C4QchAyEGwAJ

The only change is the first line :

/QGauss face_quadrature_formula(degree+2);/


2. I am using deal.II 8.4.0 now, and I think i am not using manifold 
description in my code


3. The mesh is shown as follows:

<https://lh3.googleusercontent.com/-Vs0s6So28js/WYrCbyok9mI/BSc/rfs-1D1TA8c71IQZpOEyOoU3C1idNr1DQCLcBGAs/s1600/Image.png>


Best,

Howe



在 2017年8月9日星期三 UTC+8下午3:25:20，Martin Kronbichler写道：

Dear Howe,

How did you run your simulation? From your picture, it appears
that a higher order method is worse at higher degrees than a lower
order method, which does not match with my experience. If that
were the case, nobody would use high orders. However, you need to
bring many pieces in place to really get to the benefit of the
high order method for somewhat more complicated examples such as
the flow around a cylinder. Here is a list of things to look at:

- Do you use a high-order polynomial mapping MappingQ of the same
or higher degree as the interpolation space? Do you use this
mapping in all routines that evaluate quantities, such as the
usual assembly, the computation of the lift/drag, and so on?
- Do you use a manifold description that extends into the domain?
(Look into TransfiniteInterpolationManifold.) Without, you will
not get more than third order convergence.
- Do you have a good mesh around the area of interest? Flows
around cylinders tend to be really really sensitive to the mesh
quality around the cylinder.

For the Navier-Stokes equations around the cylinder, if everything
is done right one gets significantly improved results in terms of
accuracy over the number of degrees of freedom up to degree (6,5)
(velocity,pressure). Beyond that picture is less clear. At least
with the meshes that we tried in our group it was not worth to go
beyond. You can have a look a our results in section 5.4 and Figs.
9 and 10 of this preprint:
https://arxiv.org/pdf/1706.09252.pdf
<https://arxiv.org/pdf/1706.09252.pdf>

Best,
Martin


On 09.08.2017 09:01, Howe wrote:

Dear Jaekwang

Have you solved this problem? If yes, Could pls share your
solution with us?
I am simulating a steady state flow over a cylinder, and the
drag/lift coefficient shows an unexpected trend of change as i
increase the discretization order and refine the mesh.


<https://lh3.googleusercontent.com/-Gz5932Zt6e0/WYqvbO3P-4I/BSM/6hBHn1FO5W0P7X3SPfQ4iyRaldDYzOt3QCLcBGAs/s1600/Image.png>

As is shown in the figure, the Cd increased as the cells
increased for all the discretization orders, however, for a fixed
cells, the Cd decreased as the discretization order increased.

In my opinion, to increase the order and refine the mesh should
both make the approximation more clo

Re: [deal.II] Re: Compilation fail with dealii 8.5.0 from Docker

[Martin Kronbichler]

I let someone else comment on the docker configuration and why that did
happen, but I am pretty sure that you have DEAL_II_WITH_P4EST commented
out (=disabled) in your
DEAL_II_INSTALL_DIR/include/deal.II/base/config.h which indicates that
it was eventually not enabled.

Best,
Martin


On 11.05.2017 14:47, Ashkan Dorostkar wrote:
> I don't have a dealii build directory. I tried using the docker image.
> However I see that p4est is installed.
> In the container if I do 
> env | grep DIR
> I get the following
>
> P4EST_DIR=/home/dealii/libs/p4est-1.1
>
> HDF5_DIR=/home/dealii/libs/hdf5-1.10.0-patch1
>
> DEAL_II_DIR=~/dealii-v8.5.0
>
> MUMPS_DIR=/home/dealii/libs/petsc-3.7.4
>
> METIS_DIR=/home/dealii/libs/petsc-3.7.4
>
> SUPERLU_DIST_DIR=/home/dealii/libs/petsc-3.7.4
>
> OPENCASCADE_DIR=/home/dealii/libs/oce-0.17.2
>
> TRILINOS_DIR=/home/dealii/libs/trilinos-12-8-1
>
> SLEPC_DIR=/home/dealii/libs/slepc-3.7.3
>
> SUPERLU_DIR=/home/dealii/libs/petsc-3.7.4
>
> ARPACK_DIR=/home/dealii/libs/arpack-3.4.0
>
> PARMETIS_DIR=/home/dealii/libs/petsc-3.7.4
>
> SCALAPACK_DIR=/home/dealii/libs/petsc-3.7.4
>
> PETSC_DIR=/home/dealii/libs/petsc-3.7.4
>
> OCE_DIR=/home/dealii/libs/oce-0.17.2
>
>
> On Thursday, May 11, 2017 at 2:09:01 PM UTC+2, Martin Kronbichler wrote:
>
> Dear Ashkan,
>
> the error message pops up from line 1043, which is the case when
> your deal.II configuration does not have P4EST. This is the same
> code as in 8.4.0. It seems that your deal.II configuration misses
> to find P4EST, i.e., DEAL_II_WITH_P4EST is off with the 8.5
> version. Can you check in the detailed.log file in your deal.II
> build directory whether something has changed?
>
> Best,
> Martin
>
> On 11.05.2017 12:48, Ashkan Dorostkar wrote:
>> I have attached the CMake file and the source file. This is as
>> minimal as I can make it. The complete error is as follows
>>
>> /home/dealii/app/mock/step-3.cc:In constructor ‘Test::mock::mock()’:
>>
>> /home/dealii/app/mock/step-3.cc:28:65:error: no matching function
>> for call to
>> 
>> ‘dealii::parallel::distributed::Triangulation<2>::Triangulation(ompi_communicator_t*&,
>> dealii::Triangulation<2, 2>::MeshSmoothing)’
>>
>>  Triangulation::smoothing_on_coarsening))
>>
>>  ^
>>
>> In file included from /home/dealii/app/mock/step-3.cc:1:0:
>>
>> 
>> /home/dealii/dealii-v8.5.0/include/deal.II/distributed/tria.h:1043:7:note:
>> candidate: dealii::parallel::distributed::Triangulation<dim,
>> spacedim>::Triangulation() [with int dim = 2; int spacedim = 2]
>>
>>Triangulation();
>>
>>^
>>
>> 
>> /home/dealii/dealii-v8.5.0/include/deal.II/distributed/tria.h:1043:7:note:
>>   candidate expects 0 arguments, 2 provided
>>
>> 
>> /home/dealii/dealii-v8.5.0/include/deal.II/distributed/tria.h:1037:11:note:
>> candidate:
>> dealii::parallel::distributed::Triangulation<2>::Triangulation(const
>> dealii::parallel::distributed::Triangulation<2>&)
>>
>>  class Triangulation: public
>> dealii::parallel::Triangulation<dim,spacedim>
>>
>>^
>>
>> 
>> /home/dealii/dealii-v8.5.0/include/deal.II/distributed/tria.h:1037:11:note:
>>   candidate expects 1 argument, 2 provided
>>
>> 
>> /home/dealii/dealii-v8.5.0/include/deal.II/distributed/tria.h:1037:11:note:
>> candidate:
>> 
>> dealii::parallel::distributed::Triangulation<2>::Triangulation(dealii::parallel::distributed::Triangulation<2>&&)
>>
>> 
>> /home/dealii/dealii-v8.5.0/include/deal.II/distributed/tria.h:1037:11:note:
>>   candidate expects 1 argument, 2 provided
>>
>> CMakeFiles/step-3.dir/build.make:62: recipe for target
>> 'CMakeFiles/step-3.dir/step-3.cc.o' failed
>>
>> make[2]: *** [CMakeFiles/step-3.dir/step-3.cc.o] Error 1
>>
>> CMakeFiles/Makefile2:227: recipe for target
>> 'CMakeFiles/step-3.dir/all' failed
>>
>> make[1]: *** [CMakeFiles/step-3.dir/all] Error 2
>>
>> Makefile:83: recipe for target 'all' failed
>>
>> make: *** [all] Error 2 
>>
>> On Thursday, May 11, 2017 at 11:42:18 AM UTC+2, Daniel Arndt wrote:
>>
>> Ashkan,
>>
>> I trie

Re: [deal.II] Effect of attaching SphericalManifold to all entities

[Martin Kronbichler]

Hi Praveen,

no, MappingQ does not unless you set the second boolean flag in the
constructor to enable the mapping on all cells and not just the
boundary. Even worse, when you have attached a curved manifold on the
cells and not just the boundary you will have gaps in your computational
domain. I highly recommend using MappingQGeneric instead.

Best,
Martin

P.S. We had a discussion about removing the class MappingQ and just use
MappingQGeneric always just the other day:
https://github.com/dealii/dealii/issues/3874#issuecomment-305513405 (and
subsequent comments).


On 02.06.2017 03:25, Praveen C wrote:
> If I do the following   
>
> (dim = 2)
>
> GridGenerator::hyper_shell (triangulation,
>Point(0,0),
>1.0,
>2.0);
>static const SphericalManifold boundary;
>triangulation.set_all_manifold_ids(0);
>triangulation.set_manifold (0, boundary);
>
> and use MappingQ(2), does this make all faces in the mesh to be curved 
> (except the radially oriented faces) ?
>
> 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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] AssertDimension (partition_blocks[partition], task_info.n_blocks) when calling MatrixFree::reinit() for finest level in GMG

[Martin Kronbichler]

Hi Denis,

This is a bug somewhere in the initialization of the blocks. Could you 
please send me the code to reproduce this behavior?


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Matrix-Free framework for vector problems: no assertion for components

[Martin Kronbichler]

Dear Michal,

Could you please be a bit more specific? What was the error that you 
obtained? What vector did you pass into 
FEEvaluation::read_dof_values/distribute_local_to_global? I guess that 
you handed in a simple vector, e.g. LinearAlgebra::distributed::Vector 
and got a segmentation fault. But to add the assertion I need more 
details. The issues that in principle, we do support using multiple 
components out of a scalar FE - but in that case the user must hand in a 
BlockVector or std::vector. We could catch that case in the 
machinery, though.


Thanks!

Best,
Martin


On 07.10.2017 23:00, Michał Wichrowski wrote:

Dear deal.II developers,
I've found out that in Matrix-free framework there is no assertion 
that checks if number of components of finite element matches number 
components used in cell operations. By mistake I used scalar FE 
instead vector and it took me quite a long time to figure it out. The 
program worked and even covered in reasonable number of iterations. I 
think FEEvaluation is right place to put this assertion.


Michał Wichrowski
--
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 
.

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Bug in FEEvaluationAccess::get_symmetric_gradient with Number =float

[Martin Kronbichler]

Dear Michal,

You are right, this is a bug and your fix is indeed correct. Would you 
like to create a pull request fixing it? Here are some basic 
instructions: http://dealii.org/participate.html including a video 
lecture. This way you get the proper credit for it.


Best,
Martin


On 08.10.2017 12:38, Michał Wichrowski wrote:

Dear all,
In following line:
line 4545 of /include/deal.II/matrix_free/fe_evaluation.h:
  VectorizedArray half = make_vectorized_array (0.5);

0.5 is interpreted as double and it leads to using 
make_vectorized_array. That is OK if Number=double but for any 
other type (eg. float) it results is compilation error. It should be 
replaced by:

  VectorizedArray half = make_vectorized_array (0.5);

Michał
--
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 
.

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Bug in FEEvaluationAccess::get_symmetric_gradient with Number =float

[Martin Kronbichler]

Dear Michal,

You need to create a fork of deal.II, push to your deal.II fork, and
then create the pull request from the web interface on your fork.

Best,
Martin


On 10.10.2017 13:45, Michał Wichrowski wrote:
> I've tried to create pull request, but I've got error message:
>
> mwichro@Preludio:~/lib/dealii$ git push origin 
> fix_feevaluation_float_inst
> Username for 'https://github.com': mwichro
> Password for 'https://mwic...@github.com': 
> remote: Permission to dealii/dealii.git denied to mwichro.
> fatal: unable to access 'https://github.com/dealii/dealii.git/': The
> requested URL returned error: 403
>
>
> W dniu niedziela, 8 października 2017 12:58:50 UTC+2 użytkownik Martin
> Kronbichler napisał:
>
> Dear Michal,
>
> You are right, this is a bug and your fix is indeed correct. Would
> you like to create a pull request fixing it? Here are some basic
> instructions: http://dealii.org/participate.html
> <http://dealii.org/participate.html> including a video lecture.
> This way you get the proper credit for it.
>
> Best,
> Martin
>
>
> -- 
> 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
> <mailto:dealii+unsubscr...@googlegroups.com>.
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Speeding up characteristic interpolation

[Martin Kronbichler]

Dear Praveen,

Regarding the performance of your code, there are several fundamental
problems:
- You run the constructor of Quadrature (allocates memory) and FEValues
(allocates lots of memory, evaluates a number of things you don't need)
that are both not made for use within the scenario of FEFieldFunction or
your manual loop.
- Since you're only using the shape values, you could skip all the
FEValues stuff and simply run the loop as
  solution += solution0(local_dof_indices[i]) * fe.shape_value(i, p_ref);
- At this point, the performance depends on the
transform_real_to_unit_cell (which might be slow due to Newton
iterations) and the polynomial evaluation of fe.shape_value and to some
extent the virtual function call of the latter. It depends on your
polynomial space whether this is expensive.
- You could speed up the polynomial evaluation significantly if the
polynomial basis had a way to evaluate all polynomials at once - one
could use the tensor product structure here. If you're interested, I
could give you some code to do that.

I am not aware of any generic performance optimization beyond that level
for generic mappings, but I'm also not an expert of semi-Lagrangian
methods. The works that I know of do usually work on Cartesian grids
where one can get tremendous speedups of the
transform_real_to_unit_cell. Then you spend most of the time in the
polynomial interpolation steps again - but probably a factor of 100
faster than with what you have now.

Best,
Martin


On 05.09.2017 15:18, Praveen C wrote:
> Hello Bruno
>
>> Are you sure that  FEFieldFunction.value() is the slowest part? In
>> 3D, point_inside could also be slow (you may need to do a few Newton
>> Iterations). It would be nice if you ran your code using a profiler
>> to see what is the bottleneck (building the quadrature in value(p1),
>> the evaluation itself, etc.).
>>
>
> If I comment out call to FEFieldFunction.value(), the code runs a LOT
> faster, so there is no doubt this is the slow part. I am in 2-d at
> present.
>
> I have set the correct cell by call to set_active_cell. Does the value
> function again check if the point is in this cell ?
>
> I implemented my own function to evaluate the solution as follows, but
> this is equally slow. I havent used a profiler but will timing parts
> of this below function be useful ?
>
> template
> doubleProblem::evaluate_solution
> (typenameDoFHandler::active_cell_iterator cell,
> constPoint ) const
> {
>// Transform p to unit cell
>Point p_ref = mapping.transform_real_to_unit_cell(cell, p);
>
>// Create quadrature rule
>Quadrature quadrature(p_ref);
>
>FEValues fe_values (mapping,
> fe,
> quadrature,
> update_values);
>std::vector local_dof_indices
> (fe.dofs_per_cell);
>
>fe_values.reinit(cell);
>cell->get_dof_indices (local_dof_indices);
>doublesolution = 0.0;
>for(unsignedinti=0; i   solution += solution0(local_dof_indices[i]) *
> fe_values.shape_value(i, 0);
>
>returnsolution;
> }
>
> I am processing one point at a time which could be slow. I am now
> rewriting the code so that all points falling in a cell are evaluated
> at once. 
>
> 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
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Speeding up characteristic interpolation

[Martin Kronbichler]

Hi Praveen,


- At this point, the performance depends on the 
transform_real_to_unit_cell (which might be slow due to Newton 
iterations) and the polynomial evaluation of fe.shape_value and to 
some extent the virtual function call of the latter. It depends on 
your polynomial space whether this is expensive.


Currently I am on Cartesian grids but I will use quads in future. I 
will use only FE_Q(2).


Assume you want to speed up the evaluation of FE_Q::shape_value:
If you know that you only use FE_Q(2), I would explicitly write 
down the polynomial space in a helper function and first evaluate the 1D 
polynomials in the dim coordinate directions. You could either build the 
polynomial basis through the deal.II facilities by

Polynomials::LagrangeEquidistant::generate_complete_basis(2)
https://github.com/dealii/dealii/blob/master/include/deal.II/base/polynomial.h#L345
or even hardcode the 9 coefficients of the quadratic basis in 1D 
according to

https://github.com/dealii/dealii/blob/master/source/base/polynomial.cc#L772
The latter is going to be faster because Polynomial::value has some 
overhead.


Once you have the dim arrays of 1D values, you can expand the tensor 
product. This is a series of dim nested loops where you multiply the 
respective indices. I do not know by hard if we have this pattern 
documented well somewhere, but there is some uses of this in 
FE_Q_Base::get_restriction_matrix for example. Then you also need to 
reorder according to the numbering in FE_Q. Fortunately, there is 
TensorProductPolynomials::compute,

https://github.com/dealii/dealii/blob/master/include/deal.II/base/tensor_product_polynomials.h#L115
which evaluates the whole basis. It does not contain the optimization 
mentioned above, but that would be easy to add. Since FE_Q does not 
expose the polynomial space, I'd manually construct a 
TensorProductPolynomials class and manually apply the numbering of FE_Q 
by FE_Q::get_poly_space_numbering_inverse.


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Matrix-free: Computing coefficient involving material ID

[Martin Kronbichler]

Dear Michal,


>   for(unsigned int element=0;
> elementdata->n_components_filled(cell) instead. It returns how
many of the SIMD lanes are actually filled with real data.

Best,
Martin

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Matrix-Free framework for vector problems: no assertion for components

[Martin Kronbichler]

Dear Michal,

thanks for the use case. I have now created a pull request to the
deal.II master branch that should add the appropriate assertion:

https://github.com/dealii/dealii/pull/5212

The problem was that we previously set all FE components to the single
scalar vector component, effectively duplicating the components and
duplicating the results. We now correctly identify the "single" vector
given to the evaluation routine and set nullptr for the others which
allows us to catch unintended use.

Best,
Martin


On 08.10.2017 13:19, Michał Wichrowski wrote:
>
>
> W dniu niedziela, 8 października 2017 12:55:02 UTC+2 użytkownik Martin
> Kronbichler napisał:
>
> Dear Michal,
>
> Could you please be a bit more specific? What was the error that
> you obtained? What vector did you pass into
> FEEvaluation::read_dof_values/distribute_local_to_global? I guess
> that you handed in a simple vector, e.g.
> LinearAlgebra::distributed::Vector and got a segmentation fault.
> But to add the assertion I need more details. The issues that in
> principle, we do support using multiple components out of a scalar
> FE - but in that case the user must hand in a BlockVector or
> std::vector. We could catch that case in the
> machinery, though.
>
> Thanks!
>
> Best,
> Martin
>
> Martin,
> There was no segmentation fault, program runned without any errors. I
> realised that something is wrong when I outputed results.  This is
> mine local apply,  VectorType =LinearAlgebra::distributed::Vector,
> problem (Stokes) consists of dim components of velocity and 1
> component of pressure. 
> Instead of using FESystem(FE_Q, dim) I used FE_Q for velocity.
> After fixing this it work OK, so it is just a problem with assertion.
> I'll try to make minimal code reproducing error  (or rather lack of
> error).
>   
> template 
> void
> SymGradMass<dim, degree_u,Number>
> ::local_apply (const MatrixFree<dim,Number> ,
>  VectorType  ,
>  const VectorType,
>  const std::pair _range) const
> {
>
>   typedef VectorizedArray vector_t;
>   FEEvaluation<dim,degree_u,degree_u+1+0,dim,Number> velocity (data,
> 0);  //only velocity is used here 0-th component.
>
>   for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
> {
>   velocity.reinit (cell);
>   velocity.read_dof_values (src);
> integrate_A(velocity);
>   velocity.distribute_local_to_global (dst);
>
> }
> }
>
> Michał
> -- 
> 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
> <mailto:dealii+unsubscr...@googlegroups.com>.
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Wrong result in use of FEEvaluation with different ConstraintMatrix objects

[Martin Kronbichler]

Dear Sambit,

If the result of the two cases is different and 
ConstraintMatrix::distribute() was called in both cases, I expect there 
to be some confusion regarding the indices of ghost entries. In debug 
mode, there should be a check that the parallel partitioner of the 
vector inside FEEvaluation::read_dof_values* does match with the 
expected index numbering. Did you run in debug mode? To localize the 
issue, can you check whether you called 
"matrix_free_data.initialize_dof_vector(VECTOR_NAME, 1);" in the second 
case? Note the optional argument "1" that must match with the "1" passed 
to FEEvaluation. If the issue still appears, it must be because 
ConstraintMatrix::distribute() does not do all updates. In that case, I 
would appreciate if you can give us a workable example.


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

[Martin Kronbichler]

Dear Michal,

You mean because once you apply an inhomogeneous Dirichlet condition on 
the velocity you also get a contribution because the divergence is not 
zero? You could work around that by applying a function in the whole 
domain that satisfies the Dirichlet constraints and is divergence free, 
in which case you should not get a contribution to right right hand side 
of the divergence equation but only the momentum equation. On the other 
hand, I'm not sure how easy it would be to construct such a function...


Regarding the workaround: It could work by the modification you 
suggested, but we cannot really fix that in the library because the 
operators are meant for linear operators and not affine ones. The 
problematic issue is that there are some cases where you want to apply 
an operator with a non-zero constraint and others where you don't. For 
example, all Krylov solvers expect you to provide a matrix-vector for 
the linear operator, not the affine one as you would get when you have 
inhomogeneous constraints inside the mat-vec. In particular, I do not 
really understand the part in your second post where you talk about the 
vmult_interface_{up,down} functions because those are only used for 
multigrid where the library always assumes to work on homogeneous problems.


I have not gone down this route and as far as I know, nobody else has 
previously - including the matrix-based solvers we have available in 
deal.II whose way of operation resembles the setup I described in my 
previous post, but do it inside the 
ConstraintMatrix::distribute_local_to_global call - so I cannot say how 
much work that would be. We would be happy to provide a solution that is 
generic if you find one, but I'm not sure I will be able to help much 
along this direction.


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Transfinite Interpolation

[Martin Kronbichler]

Hi Narenda and David,

Just to add on top of what David Wells said, Luca Heltai and I once had
the idea of showing this manifold in a new tutorial program (where we
combine it with MappingFEField to show the performance impact). Now
step-60 does contain some MappingFEField, but on a different topic, so I
think we should simply create such a program soon. I will open an issue
for that.

Best,
Martin


On 23.05.2018 23:32, David Wells wrote:
> Hi Narendra
>
> We should probably put something like this in step-49. Here is a bit
> of code that shows how to use the TFI manifold class:
>
> |
> #include
> #include
> #include
> #include
>
> intmain()
> {
>   usingnamespacedealii;
>
>   Triangulation<2>tria;
>   GridGenerator::hyper_cube_with_cylindrical_hole (tria,0.05,0.1);
>   consttypes::manifold_id tfi_id =1;
>   for(constautocell :tria.active_cell_iterators())
>     {
>       // set all non-boundary manifold ids to the new TFI manifold id.
>       cell->set_manifold_id(tfi_id);
>       for(unsignedintface_n =0;face_n ::faces_per_cell;
>            ++face_n)
>         {
>           if(!cell->face(face_n)->at_boundary())
>             cell->face(face_n)->set_manifold_id(tfi_id);
>         }
>     }
>   TransfiniteInterpolationManifold<2>inner_manifold;
>   inner_manifold.initialize(tria);
>   tria.set_manifold(1,inner_manifold);
>
>   tria.refine_global(2);
>
>   GridOutgo;
>   std::ofstream out("grid.eps");
>   go.write_eps(tria,out);
> }
> |
>
>
> This creates a nice blend between the inner circle and the outer
> square. From here you can do something like what is done in step-49 to
> generate the rest of the grid from Cartesian pieces.
>
> If you are using 8.5 instead of 9.0 you will need to make all Manifold
> classes (including putting your own PolarManifold on the circle)
> before the Triangulation.
>
> I hope this helps!
>
> Thanks,
> David Wells
>
>
>
> On Monday, May 21, 2018 at 7:06:19 PM UTC-4, Narendra Nanal wrote:
>
> Hello All,
>
> I am trying to run benchmark case of 2D flow around a cylinder.
> Here are the details- 
> 
> http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark1_re20.html
> 
> .
> My question is how to use Transfinite Interpolation class to
> generate a mesh with a smooth transition. Thank you. 
>
> Regards,
> Narendra 
>
> -- 
> 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
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Add image to step

[Martin Kronbichler]

Hi Giovanni,

I do have access to the svn for the homepage, so you can send them to me
if you want (and nobody else was quicker in replying).

Best,
Martin


On 02.05.2018 16:22, giovannialze...@hotmail.it wrote:
> Hi everyone!
>    I am working with Luca on step 60, it would be nice to add a couple
> of images. Looking at the other examples the images are uploaded on
> the deal.ii website, and embedded using things like:
>
> img
> src="https://www.dealii.org/images/steps/developer/step-1.grid-1.png;
> alt=""
>
>
> Temporarily I shall put them on my sissa website, so that I can easily
> do the next push but, who can I send the images to? Thank you!
>
> Best,
> Giovanni
> -- 
> 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
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Compute diagonal with Matrix-Free on adaptive meshes

[Martin Kronbichler]

Dear Daniel,

the problem is, as far as I can tell, the fact that once you assemble 
into a matrix and once into a vector. The paper


K. Kormann: A Time-Space Adaptive Method for the Schrödinger Equation, 
Commun. Comput. Phys. 2016, doi: 10.4208/cicp.101214.021015a


describes this effect in section 5.3. Can you check there if this is 
what you see?


Best,
Martin


On 02.05.2018 20:17, Daniel Jodlbauer wrote:
But thats exactly my point, the error occurs in the dofs which 
constrain a hanging node, not the hanging nodes (dofs 5 and 7) itself. 
I agree that the constrained dofs 5 and 7 can have arbitrary values.
I will check whether the solution is going to be different in any 
case. I was just afraid that wrong values on the diagonal could cause 
problems in the MGSmoother afterwards.


On Wednesday, May 2, 2018 at 7:38:29 PM UTC+2, Wolfgang Bangerth wrote:


Daniel,

> To verify my MatrixFree implementation, I compared its
application to
> the classical matrix-vector multiplication (call it matrix *A*).
> This is done by computing the matrix *M* of the operator *MF* by
> applying it to all unit-vectors.
>
> However, when I compute the diagonal in the same way as
LaplaceOperator
> does it (copied it), I get values different from the assembled
diagonal
> once I have hanging nodes.

The rows and columns corresponding to hanging nodes are empty with
the
exception of the diagonal entry -- for which we use a value that
has the
correct order of magnitude, but that is otherwise unspecified. In
other
words, what diagonal value you have there is unimportant as long
as it
is nonzero because these degrees of freedom don't couple with all
of the
other degrees of freedom, and as long as you overwrite the values
of the
computed solution through ConstraintMatrix::distribute().

It is not surprising to me that you get different diagonal entries
when
you use two different methods. The question is whether you get a
different *solution* vector after distributing to hanging nodes.

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 
.

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Petrov-Galerkin stabilization approach to solve Navier-Stoke equations

[Martin Kronbichler]

Hi Jie,

>   Thank you for your hints on my previous question. Now I have
> extended step-57 to be time-dependent. I added a u,t term which is
> approximated using backward Euler, and all the remaining terms are
> treated implicitly. Newton iteration is still used at every time step.
> The system becomes
>   (A + M/dt)U = F - (\delta v,  (u^{k+1} - u^k)/dt)_\Omega, BU = P,
> which is similar to what we have in step-57 except the LHS(0, 0) and
> RHS(0) are changed. 
>
>   I'm still using the same preconditioner developed in step-57. So far
> it works the simple case of lid-driven cavity flow. My questions are: 
>   1. It doesn't make much sense to keep using the preconditioner
> developed for steady NSE here, are there any simple ways to adapt it
> for time-dependent problems?

You should check the literature. In case your time step is not too
large, the mass matrix will give a significant contribution to the
system matrix. This is good because one knows how to approximate the
Schur complement for a velocity mass matrix - it is simply a pressure
Poisson matrix, with Dirichlet conditions on that operator where the
velocity has Neumann and Neumann conditions where you specify the flow.
There are several variants of preconditioners in the literature. One is
Cahouet-Chabard 1988 that simply uses the inverse Poisson matrix and the
inverse mass matrix to take the mass and viscous contributions in
velocity into account. I have used it in a recent publication:
https://doi.org/10.1177/1094342016671790 - check also literature cited
there, this mailing list is not the right place to ask for a literature
review.

Then, there are also pressure reaction-convection-diffusion operators
that also take convection into account. I have not seen them a lot for
the full block system in the time-dependent case, probably because the
range where they are really better than Cahouet-Chabard yet the Schur
complement approximation stays good enough is quite narrow - especially
once you also want to have good yet cheap approximation of the velocity
matrix.

>   2. Instead of working on better preconditioners, can I simply treat
> the convection term explicitly (using the value at last time step)?
> The reason is that doing so would make the LHS of the equation
> symmetric  (I think), which should be easy to solve. This seems to be
> much easier.
Yes, this is Cahouet-Chabard and indeed very efficiently done. But you
get a CFL limit which may or may not be limiting for you.

Best,
Martin

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Calculating lift and drag

[Martin Kronbichler]

Hi Jie,

One needs to move some of the vertices around. I have some code that
does generate this geometry in one of my projects:

https://github.com/kronbichler/adaflo/blob/master/tests/flow_past_cylinder.cc

(check the create_triangulation function)

The mesh that comes out from this code is accurate but not optimal yet
it is not fine enough around the cylinder as compared to be bulk of the
domain, at least for some of the benchmark test cases. You can find a
picture which gives better convergence results in our recent publication
(Fig 11):

http://www.sciencedirect.com/science/article/pii/S0021999117306915

However, I do not have that code online right now and it takes too much
time to look it up. But it should be easy enough to compose it of
several triangulation, adding the transfinite interpolation manifold to
the transition from polar->straight elements that we recently introduced
in the deal.II and that is mentioned in that paper.

Best,
Martin


On 20.10.2017 01:15, Jie Cheng wrote:
> Hi Lucas
>
>
> I am trying to replicate the results here
> 
> 
> of the benchmark problem of flow around a cylinder (circle in 2D).
> Here is my code so far:
>
>
>   It's been a long time since you posted this question. But if you are
> still reading the threads, could you please tell me how you managed to
> generate a non-symmetric triangulation in the benchmark case? I tried
> to use one hypercube_with_cylindrical_hole and two hyper_rectangles to
> construct it, but the result looks ugly upon global refinement.
>
>Thank you.
>
> Jie
> -- 
> 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
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Calculating lift and drag

[Martin Kronbichler]

Hi Jie,

One needs to move some of the vertices around. I have some code that
does generate this geometry in one of my projects:

https://github.com/kronbichler/adaflo/blob/master/tests/flow_past_cylinder.cc

(check the create_triangulation function)

The mesh that comes out from this code is accurate but not optimal yet
it is not fine enough around the cylinder as compared to be bulk of the
domain, at least for some of the benchmark test cases. You can find a
picture which gives better convergence results in our recent publication
(Fig 11):

http://www.sciencedirect.com/science/article/pii/S0021999117306915

However, I do not have that code online right now and it takes too much
time to look it up. But it should be easy enough to compose it of
several triangulation, adding the transfinite interpolation manifold to
the transition from polar->straight elements that we recently introduced
in the deal.II and that is mentioned in that paper.

Best,
Martin


On 20.10.2017 01:15, Jie Cheng wrote:
> Hi Lucas
>
>
> I am trying to replicate the results here
> 
> 
> of the benchmark problem of flow around a cylinder (circle in 2D).
> Here is my code so far:
>
>
>   It's been a long time since you posted this question. But if you are
> still reading the threads, could you please tell me how you managed to
> generate a non-symmetric triangulation in the benchmark case? I tried
> to use one hypercube_with_cylindrical_hole and two hyper_rectangles to
> construct it, but the result looks ugly upon global refinement.
>
>Thank you.
>
> Jie
> -- 
> 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
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

[Martin Kronbichler]

Dear Michal,

This is expected: the matrix-free operator evaluation cannot apply
non-homogeneous boundary conditions while solving (at least I have never
figured out how to do that). To solve such a problem, you need to bring
the non-homogeneous part on the right hand side first. I often solve
this by first creating a vector that is filled with the Dirichlet values
and apply the operator on that (without setting the zero constraints)
and then solve for an increment that has zero boundary conditions. I
attach an example which is an extension of step-37 (and uses
coefficients described here, if you want some more background:
http://www.it.uu.se/research/publications/reports/2017-006/2017-006-nc.pdf

In the program, you will find that LaplaceOperator not only contains a
vmult() function, but also a compute_residual() function that reads the
data in from an FEEvaluation object without Dirichlet conditions, due
the operation, and assembles the right hand side into an FEEvaluation
object with Dirichlet conditions (-> local_compute_residual()).

Please let me know in case you have questions.

Best,
Martin

On 26.10.2017 22:07, Michał Wichrowski wrote:
> Dear deal.II developers,
> I think I found problem with non-homogeneus boundary conditions.  
> By changing:
> VectorTools::interpolate_boundary_values (dof_handler,
> 0,
> Functions::ZeroFunction(),
> constraints);
> to
> VectorTools::interpolate_boundary_values (dof_handler,
> 0,
> Functions::ConstantFunction(5),
> constraints);
>
>
> The results shown on included picture were outputted. I've tried to
> figure out what is going on, and it looks like matrix-free framework
> is solving problem with homogeneous constrains and then Dirichlet
> nodes are overwritten.  
> -- 
> 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
> <mailto:dealii+unsubscr...@googlegroups.com>.
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 2009 - 2017 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -----

 *
 * Authors: Katharina Kormann, Martin Kronbichler, Uppsala University,
 * 2009-2012, updated to MPI version with parallel vectors in 2016
 */


#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 


namespace Step37
{
  using namespace dealii;


  const unsigned int degree_finite_element = 3;
  const unsigned int dimension = 3;

  typedef double number;
  typedef float level_number;


  // ==
  // Functions
  // ==

  template 
  class Solution : public Function
  {
  public:
Solution ()  : Function() {}

virtual double value (const Point   ,
  const unsigned int  component = 0) const;

virtual Tensor<1,dim> gradient (const Point   ,
const unsigned int  component = 0) const;


virtual Tensor<1,dim,VectorizedArray> gradient (const 
Point<dim,VectorizedArray>  ,
const unsigned int  
component = 0) const;

virtual double laplacian (const Point   ,
  const unsigned int  component = 0) const;

virtual VectorizedAr

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

[Martin Kronbichler]

Hi Bruno,

I definitely agree - we should mention this in step-37. I will open an
issue since I will not get to do it today.

Best,
Martin


On 27.10.2017 14:44, Bruno Turcksin wrote:
> Martin,
>
> Could you add this explanation to the documentation and/or to a
> tutorial. I think mentioning it step-37 would be great.
>
> Best,
>
> Bruno
>
> On Friday, October 27, 2017 at 3:57:09 AM UTC-4, Martin Kronbichler
> wrote:
>
> Dear Michal,
>
> This is expected: the matrix-free operator evaluation cannot apply
> non-homogeneous boundary conditions while solving (at least I have
> never figured out how to do that). To solve such a problem, you
> need to bring the non-homogeneous part on the right hand side
> first. I often solve this by first creating a vector that is
> filled with the Dirichlet values and apply the operator on that
> (without setting the zero constraints) and then solve for an
> increment that has zero boundary conditions. I attach an example
> which is an extension of step-37 (and uses coefficients described
> here, if you want some more background:
> http://www.it.uu.se/research/publications/reports/2017-006/2017-006-nc.pdf
> 
> <http://www.it.uu.se/research/publications/reports/2017-006/2017-006-nc.pdf>
>
> In the program, you will find that LaplaceOperator not only
> contains a vmult() function, but also a compute_residual()
> function that reads the data in from an FEEvaluation object
> without Dirichlet conditions, due the operation, and assembles the
> right hand side into an FEEvaluation object with Dirichlet
> conditions (-> local_compute_residual()).
>
> Please let me know in case you have questions.
>
> Best,
> Martin
>
> On 26.10.2017 22:07, Michał Wichrowski wrote:
>> Dear deal.II developers,
>> I think I found problem with non-homogeneus boundary conditions.  
>> By changing:
>> VectorTools::interpolate_boundary_values (dof_handler,
>> 0,
>> Functions::ZeroFunction(),
>> constraints);
>> to
>> VectorTools::interpolate_boundary_values (dof_handler,
>> 0,
>> Functions::ConstantFunction(5),
>> constraints);
>> The results shown on included picture were outputted. I've tried
>> to figure out what is going on, and it looks like matrix-free
>> framework is solving problem with homogeneous constrains and then
>> Dirichlet nodes are overwritten.  
>> -- 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
>> <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
>> <mailto:dealii+unsubscr...@googlegroups.com>. For more options,
>> visit https://groups.google.com/d/optout
>> <https://groups.google.com/d/optout>. 
>
> -- 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
> <mailto:dealii+unsubscr...@googlegroups.com>. For more options, visit
> https://groups.google.com/d/optout. 

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Possible bug in the Jacobi preconditioner for MatrixFreeOperators::LaplaceOperator for vector fields

[Martin Kronbichler]

Dear Stephen,

You are absolutely right, the value of dofs_per_cell is simply wrong in 
the vector-valued case. I have been hesitant to fix it because there are 
some downstream projects using it (mostly mine, though), but I guess it 
is better to switch to the correct notation now rather than causing more 
trouble. You can use your fix right now - I will open an issue and fix 
it (probably next week).


Best,
Martin


On 27.10.2017 21:58, Stephen DeWitt wrote:

Hello,
I've been trying to refactor my code to use the new 
MatrixFreeOperators, but I've run into a problem trying to use the 
Jacobi preconditioner with the MatrixFreeOperators::LaplaceOperator 
with a vector-valued variable.


In short, I wrote a short code to solve a simple 2D Poisson problem. 
For a scalar field, it works well both when using the identity matrix 
preconditioner and the Jacobi preconditioner. However, for a vector 
field, it works when using the identity matrix preconditioner but 
won't converge with the Jacobi preconditioner.


Aside from all of the standard book-keeping (creating an FESystem with 
multiple components, applying vector-valued ICs and BCs, etc.), my 
only change between the scalar case and the vector case is the 
template parameter for the MatrixFreeOperators::LaplaceOperator object.


Digging into it, I found that approximately every other element of the 
diagonal vector was zero, excluding a few that I believe are from the 
BCs (that is, in "compute_diagonal()", half of the elements in 
"inverse_diagonal_vector" are zero after the cell loop but before the 
reciprocal is taken). In the scalar case, all of the diagonal elements 
are non-zero, as one would expect. The zero elements of the diagonal 
seem to come from the "local_diagonal_cell" method, where 
"phi.dofs_per_cell" and "phi.tensor_dofs_per_cell" don't change 
between the scalar and vector case (4 in this case, with linear 
elements), when I'd expect the dofs per cell in the vector case to be 
the number of components times the dofs per cell in the scalar case.


If I multiply  "phi.dofs_per_cell" and "phi.tensor_dofs_per_cell" by 
the number of components everywhere in the "local_diagonal_cell" 
method, the preconditioner works (see the highlighted edits at the 
bottom of the post).


So, does anyone have an idea of what's going on here? Is this a bug, 
or is there an extra step I'm missing to use MatrixFreeOperators with 
a vector that I'm compensating for?


Thanks,
Steve

|
template
void
CustomLaplaceOperator::
  local_diagonal_cell 
(constMatrixFree::value_type>,

VectorType,
constunsignedint&,
conststd::pair_range)const
{
typedeftypenamedealii::MatrixFreeOperators::Base::value_type 
Number;
FEEvaluationphi 
(data,this->selected_rows[0]);


for(unsignedintcell=cell_range.first;cell.

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: order of cells in two MatrixFree objects on the same Triangulation and Quadrature but different FE spaces

[Martin Kronbichler]

Dear Denis,


> p.s. the mapping is the same for both cases.
>
> On Monday, January 29, 2018 at 4:16:48 PM UTC+1, Denis Davydov wrote:
>
> Dear all,
>
> I am playing around with something in the following setup:
> 1. I have a single p::d::Tria
> 2. two different DoFHandler's using scalar FE_Q of different
> degrees
> 3. The same quadrature rule (QGauss) used for both cases
> 4. Two different MatrixFree objects
>
> I am not 100% sure, but it looks like the iteration over SIMD
> chunks of cells via two different matrix-free
> objects actually leads to the different sequence of cells. I see
> that when I try to calculate some integral quantities
> of my quadrature point data.
>

I guess that's possible. Once there are different DoFHandler objects,
the algorithm that collects cells and puts them into batches might
change. The documentation does state the following (in FEEvaluation):
https://www.dealii.org/developer/doxygen/deal.II/classFEEvaluation.html
section on "Handling several integration tasks and data storage in
quadrature points" towards the end of the general class documentation.
It says "For the combination of different fields, including different
solution vectors that come from different time steps, it is mandatory
that all FEEvaluation

objects share the same MatrixFree

object."

Did you expect anything in addition?

Best,
Martin

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Vertex duplication in DataOut for discontinuous fields - reference to that?

[Martin Kronbichler]

Dear users and developers,

I have a question regarding our output strategy: When we output our 
solution fields, we always create separate vertices for each element 
because this allows us to output discontinuous fields over element edges, 
for example. We also have an FAQ entry about this: 
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#the-graphical-output-files-dont-make-sense-to-methey-seem-to-have-too-many-degrees-of-freedom

What I would like to know is if anyone knows if that scheme was ever 
written down somewhere. I'm looking for a reference I can cite. I looked 
through a few documents, including the 2007 deal.II paper, but I couldn't 
find anything. I have also seen some other FE packages doing the same, but 
so I'd also take those references, if we don't have our own. (Or I simply 
link to our FAQ, even though I prefer to give people credit through a 
Scopus reference if available.)

Best,
Martin

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] New deal.II 8.5.1 is 20% slower than deal.II 8.0.0

[Martin Kronbichler]

In general, we strive to make deal.II faster with new releases, and for 
many cases that is also true as I can confirm from my applications. I 
have ran step-23 on release 8.0 as well as the current development 
sources and I can confirm that the new version is slower on my machine. 
If I disable output of step-23, I get a run time of 4.7 seconds for 
version 8.0 and 5.3 seconds for the current version. After some 
investigations I found out that while some solver-related operations got 
faster indeed (the problem with 16k dofs is small enough to run from L3 
cache in my case), we are slower in the FEValues::reinit() calls. This 
call appears in VectorTools::create_right_hand_side() and the 
VectorTools::interpolate_boundary_values in the time loop. The reason 
for this is that we nowadays call 
"MappingQGeneric::compute_mapping_support_points" also for the bilinear 
mapping MappingQ1, which allocates and de-allocates a vector. While this 
is uncritical on higher order mappings, in 2D with linear shape 
functions the time spent there is indeed not negligible. This is indeed 
unfortunate for your use case, but I want to stress that the changes 
were made in the hope to make that part of the code more reliable. 
Furthermore, those parts of the code are not performance critical and 
not accurately tracked. It is a rather isolated issue that got worse 
here, so from this single example one definitely not say that we are 
going the wrong direction as a project.


While there are plenty of things I could imagine to make this particular 
case more efficient in the application code, way beyond the performance 
of what the version 8.0 provided - note that I would not write the code 
like that if it were performance critical - the only obvious thing is 
that we could try to work around the memory allocations by not returning 
a vector in MappingQGeneric::compute_mapping_support_points but rather 
fill an existing array in 
MappingQGeneric::InternalData::mapping_support_points. Nobody of us 
developers has this high on the priority list right now, but we would 
definitely appreciate if some of our users, like you, wants to look into 
that. I could guide you to the right spots.


Best regards,
Martin


On 26.12.2017 21:22, drgulev...@gmail.com wrote:

Yes, the two are attached. The key lines from their diff result:

$ diff detailed.log-v8.1.0 detailed.log-v8.5.1
...
< #  Compiler flags used for this build:
< #    CMAKE_CXX_FLAGS:  -pedantic -fpic -Wall 
-Wpointer-arith -Wwrite-strings -Wsynth -Wsign-compare -Wswitch 
-Wno-unused-local-typedefs -Wno-long-long -Wno-deprecated 
-Wno-deprecated-declarations -std=c++11 -Wno-parentheses -Wno-long-long
< #    DEAL_II_CXX_FLAGS_RELEASE:    -O2 -funroll-loops 
-funroll-all-loops -fstrict-aliasing -Wno-unused

---
> #  Base configuration (prior to feature configuration):
> #    DEAL_II_CXX_FLAGS:    -pedantic -fPIC -Wall -Wextra 
-Wpointer-arith -Wwrite-strings -Wsynth -Wsign-compare -Wswitch 
-Woverloaded-virtual -Wno-long-long -Wno-deprecated-declarations 
-Wno-literal-suffix -std=c++11
> #    DEAL_II_CXX_FLAGS_RELEASE:    -O2 -funroll-loops 
-funroll-all-loops -fstrict-aliasing -Wno-unused-local-typedefs

18c19
< #    DEAL_II_LINKER_FLAGS: -Wl,--as-needed -rdynamic 
-pthread

---
> #    DEAL_II_LINKER_FLAGS: -Wl,--as-needed -rdynamic 
-fuse-ld=gold

...
> #    BOOST_CXX_FLAGS = -Wno-unused-local-typedefs
...
> #  ( DEAL_II_WITH_BZIP2 = OFF )
> #    DEAL_II_WITH_CXX11 = ON
> #  ( DEAL_II_WITH_CXX14 = OFF )
> #  ( DEAL_II_WITH_GSL = OFF )
...
> #    THREADS_CXX_FLAGS = -Wno-parentheses
> #    THREADS_LINKER_FLAGS = -pthread


On Tuesday, December 26, 2017 at 10:10:44 PM UTC+3, Matthias Maier wrote:

Would you mind sending us the "detailed.log" files?

Best,
Matthias


On Tue, Dec 26, 2017, at 12:35 CST, drgul...@gmail.com
 wrote:

> Thanks. This is strange as I still get 15-20% consistently
better results
> in favor of older versions
>  on three different machines already. Two more studies on other
systems
> attached below.
>
> TEST: Step-23 (integration time modified from 5 to 150, output
suppressed)
> CMAKE_BUILD_TYPE: "Release".
>
> MACHINE 1: Intel(R) Core(TM) i7-6700K CPU @ 4.00GHz
> $ cat /etc/redhat-release
> CentOS Linux release 7.4.1708 (Core)
> $ gcc --version
> gcc (GCC) 4.8.5 20150623 (Red Hat 4.8.5-16)
>
> deal v8.1.0 (built and installed from source):
>
> $ time ./step-23
> real    1m23.768s
> user    5m46.080s
> sys    0m4.079s
>
> deal v8.5.1 (built and installed from source):
>
> $ time ./step-23
> real    1m42.416s
> user    5m37.018s
> sys    0m4.340s
>
> MACHINE 2: Intel(R) Xeon(R) CPU X5690  @ 3.47GHz
> $ lsb_release -a
> Description:    Ubuntu 14.04.5 LTS
> $ gcc --version
> gcc (Ubuntu 

Re: [deal.II] New deal.II 8.5.1 is 20% slower than deal.II 8.0.0

[Martin Kronbichler]

Dear Dmitry,

Thanks for the info. This sounds interesting. May I ask some more 
details: Are you using Q1 elements in your 2D model (this is the element 
used in step-23)? I have not looked into your history model, but I guess 
it can be adapted. Since your code is limited by 
FEValues::get_function_values, I suggest you take a look at the step-48 
tutorial program. That code does indeed have performance-tuned 
components. It requires a somewhat different approach to write the 
loops, but it should pay off for time-dependent cases, I expect a 
speedup of a factor 4-6 in 2D for linear elements for the part 
corresponding to FEValues::get_function_values and FEValues::reinit. I 
would be happy to guide you more closely if this sounds applicable.


Best,
Martin


On 27.12.2017 14:16, drgulev...@gmail.com wrote:

Thank you!

Some guidance how I could optimize the code would be appreciated. I am 
using deal.II for solving a time-dependent nonlinear 2D problem (sort 
of sine-Gordon, but a more advanced model which includes a history 
dependence, https://github.com/drgulevich/mitmojco). Most of the time 
the deal.II code spends in:


1. fe_values.get_function_values -- most of the wall time (70%)
2. fe_values.reinit -- less often
3. CG solver -- even less often

Kind regards,
Dmitry

On Wednesday, December 27, 2017 at 11:49:42 AM UTC+3, Martin 
Kronbichler wrote:


In general, we strive to make deal.II faster with new releases,
and for many cases that is also true as I can confirm from my
applications. I have ran step-23 on release 8.0 as well as the
current development sources and I can confirm that the new version
is slower on my machine. If I disable output of step-23, I get a
run time of 4.7 seconds for version 8.0 and 5.3 seconds for the
current version. After some investigations I found out that while
some solver-related operations got faster indeed (the problem with
16k dofs is small enough to run from L3 cache in my case), we are
slower in the FEValues::reinit() calls. This call appears in
VectorTools::create_right_hand_side() and the
VectorTools::interpolate_boundary_values in the time loop. The
reason for this is that we nowadays call
"MappingQGeneric::compute_mapping_support_points" also for the
bilinear mapping MappingQ1, which allocates and de-allocates a
vector. While this is uncritical on higher order mappings, in 2D
with linear shape functions the time spent there is indeed not
negligible. This is indeed unfortunate for your use case, but I
want to stress that the changes were made in the hope to make that
part of the code more reliable. Furthermore, those parts of the
code are not performance critical and not accurately tracked. It
is a rather isolated issue that got worse here, so from this
single example one definitely not say that we are going the wrong
direction as a project.

While there are plenty of things I could imagine to make this
particular case more efficient in the application code, way beyond
the performance of what the version 8.0 provided - note that I
would not write the code like that if it were performance critical
- the only obvious thing is that we could try to work around the
memory allocations by not returning a vector in
MappingQGeneric::compute_mapping_support_points but rather fill an
existing array in
MappingQGeneric::InternalData::mapping_support_points. Nobody of
us developers has this high on the priority list right now, but we
would definitely appreciate if some of our users, like you, wants
to look into that. I could guide you to the right spots.

Best regards,
Martin


On 26.12.2017 21:22, drgul...@gmail.com  wrote:

Yes, the two are attached. The key lines from their diff result:

$ diff detailed.log-v8.1.0 detailed.log-v8.5.1
...
< #  Compiler flags used for this build:
< #    CMAKE_CXX_FLAGS:  -pedantic -fpic -Wall
-Wpointer-arith -Wwrite-strings -Wsynth -Wsign-compare -Wswitch
-Wno-unused-local-typedefs -Wno-long-long -Wno-deprecated
-Wno-deprecated-declarations -std=c++11 -Wno-parentheses
-Wno-long-long
< #    DEAL_II_CXX_FLAGS_RELEASE:    -O2 -funroll-loops
-funroll-all-loops -fstrict-aliasing -Wno-unused
---
> #  Base configuration (prior to feature configuration):
> #    DEAL_II_CXX_FLAGS:    -pedantic -fPIC -Wall
-Wextra -Wpointer-arith -Wwrite-strings -Wsynth -Wsign-compare
-Wswitch -Woverloaded-virtual -Wno-long-long
-Wno-deprecated-declarations -Wno-literal-suffix -std=c++11
> #    DEAL_II_CXX_FLAGS_RELEASE:    -O2 -funroll-loops
-funroll-all-loops -fstrict-aliasing -Wno-unused-local-typedefs
18c19
< #    DEAL_II_LINKER_FLAGS: -Wl,--as-needed -rdynamic -pthread
---
> #    DEAL_II_LINKER_FLAGS: -Wl,--as-

Re: [deal.II] Re: Adding an independent DoF to a Trilinos system

[Martin Kronbichler]

I did some digging into the Trilinos reinit_matrix, and I've found 
that all of my lost time is coming from the operation


   graph->FillComplete(input_col_map, input_row_map)

This step is accounting for as much as 20 seconds per call for a 
1x200,000 matrix on 4 processors. I also found that this bottleneck 
is non-existent running in serial. There is a check earlier in 
reinit_matrix which calls graph.reset differently for 1 processor vs. 
multiple processors. The comments suggest that there's an explicit 
reason for treating the mpi case this way, but that treatment seems 
to be causing the problem for my case.

I think Martin Kronbichler might have to tell you why this is.

But if you're adventurous, take a look how the FillComplete function 
in Trilinos looks like and whether you could possibly find the same 
kinds of savings there that I found in deal.II. I suspect that it is 
this function:


https://github.com/trilinos/Trilinos/blob/master/packages/epetra/src/Epetra_CrsGraph.cpp#L974 

I'm not sure where exactly the time is spent, but I could imagine that 
the most likely spot is where the code tries to sort an array of 200,000 
entries with shell_sort which, if I remember it correctly, does not have 
a short way out if the array is already sorted and is of O(n^2) 
complexity. Not sure how one would approach this - one could suggest a 
patch, but it might take a little while until the Trilinos folks react 
(Epetra is not as actively developed any more, so there is little 
resources for it).


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Status of matrix-free CUDA support

[Martin Kronbichler]

Steve,

Just to add on what Bruno said: Explicit time integration typically 
requires only a subset of what you need for running an iterative solver 
(no decision making like in convergence test of conjugate gradient, no 
preconditioner), so running that should be pretty straight-forward.


Best,
Martin

On 05.11.18 21:06, Bruno Turcksin wrote:

Steve

On Monday, October 29, 2018 at 10:43:48 AM UTC-4, Stephen DeWitt wrote:

So far I've read through the "CUDA Support" issue
, the "Roadmap for
inclusion of GPU implementation of matrix free in Deal.II" issue
, the Doxygen
documentation for classes in the CUDAWrappers namespace
,
andthe manuscript by Karl Ljungkvist
.
Are there any other pages I should be looking at?

No that's pretty much it.

My understanding from these pages is that deal.II has partial
support for using CUDA with matrix free calculations. Currently,
calculations can be done with scalar variables (but not vector
variables) and adaptive meshes.

That's correct with the weird exception that you cannot impose 
constraints (and thus have hanging nodes) if dim equals 2 and 
fe_degree is even. There is a bug in that case but we don't know why.


A few (somewhat inter-related) questions:
1). Do all of the tools exist to create a GPU version of step-48?
Has anyone done so?

I think so but nobody as tried so I am not 100% sure. Note that we 
don't have MPI support yet.


2). What exactly would be involved in creating a GPU version of
step-48? Is it just changing the CPU Vector, MatrixFree, and
FEEvaluation classes to their GPU counterparts, plus packaging
some data (plus local apply functions?) into a
CUDAWrappers::MatrixFree< dim, Number >::Data


 struct?

Yes, pretty much. You can take a look at the matrix_free_matrix_vector 
tests here . 
These tests are doing a matrix-vector multiplication using an 
Helmholtz operator. As you can see here 
 
the main difference between the GPU and the CPU code is that the body 
loop over the quadrature points needs to be encapsulated in a functor 
when using GPU. The tests are based on the CPU tests matrix_vector 
here . 
So you should have a good idea of the modifications required to use 
the GPU code.


3). Most of the discussions seemed to revolve around linear
solves. For something like step-48 with explicit updates, will the
current paradigm work well? Or would that require shuttling data
between the GPU and CPU every time step, causing too much
overhead? (I know that in general GPUs can work very well for
explicit codes.)

That should work fine. Most people are interested in using matrix-free 
with a solver but I don't see a reason why it would be slow in step 
48. The communication between the CPU and the GPU should be limited in 
step-48.


Best,

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

For more options, visit https://groups.google.com/d/optout.


--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Linking error while compiling examples: "undefined reference to `dealii::VectorizedArray::n_array_elements'"

[Martin Kronbichler]

Dear Maxi,

This variable is populated by a series of checks on your particular
system. It happens when the configuration stage prints the following
lines, and I assume that this is what you see:

-- Performing Test DEAL_II_HAVE_SSE2
-- Performing Test DEAL_II_HAVE_SSE2 - Failed
-- Performing Test DEAL_II_HAVE_AVX
-- Performing Test DEAL_II_HAVE_AVX - Failed
-- Performing Test DEAL_II_HAVE_AVX512
-- Performing Test DEAL_II_HAVE_AVX512 - Failed
-- Performing Test DEAL_II_HAVE_ALTIVEC
-- Performing Test DEAL_II_HAVE_ALTIVEC - Failed
-- Performing Test DEAL_II_HAVE_OPENMP_SIMD
-- Performing Test DEAL_II_HAVE_OPENMP_SIMD - Success

Since you are on x86-64, you should at least have a level of 1 (anyway,
also level 0 should not result in a linker error and we need to
investigate it).

Can you please post the file CMakeFiles/CMakeError.log which contains
the output of all checks that ended up being "failed" (i.e., the checks
our configuration system is doing to determine a stable configuration
that does not fail)? This way, I can look into why that failed on your
system and it might give me some hint...

Best,
Martin


On 10.10.2018 08:24, 'Maxi Miller' via deal.II User Group wrote:
> Yes, I think that is the reason for that behavior. What can I do to
> set the level to >= 1?
> Thanks!
>
> Am Mittwoch, 10. Oktober 2018 01:14:06 UTC+2 schrieb David Wells:
>
> Hey there,
>
> I think this may be actually a bug in how we recently fixed a
> related problem:
>
> https://github.com/dealii/dealii/pull/7198
> 
>
> Could you please post the detailed.log file from when you
> configured the library? I think the value for
> DEAL_II_COMPILER_VECTORIZATION_LEVEL may be wrong (or, at least,
> it is causing problems somewhere else).
>
> Thanks,
> David Wells
>
> On Tuesday, October 9, 2018 at 6:29:25 AM UTC-4, Maxi Miller wrote:
>
> I tried to compile deal.II on a new computer, but I got
> multiple linking errors for some of the examples:
> |
> 
> /usr/lib64/gcc/x86_64-suse-linux/8/../../../../x86_64-suse-linux/bin/ld:../lib/libdeal_II.g.so.9.1.0-pre:undefinedreference
> to `dealii::VectorizedArray::n_array_elements'
> 
> /usr/lib64/gcc/x86_64-suse-linux/8/../../../../x86_64-suse-linux/bin/ld:
> ../lib/libdeal_II.g.so.9.1.0-pre: undefined reference to
> `dealii::VectorizedArray::n_array_elements'
> |
>
> The examples in question are all debug versions of the
> examples, with example number 2, 3, 24, 29, 36, 44, 45, 48,
> 50, 52, 53, 55, 56, 59, 60. I used the following external
> libraries/options:
> ADOL-C, Boost, C++14/17, GSL, HDF5, Lapack, Metis, MPI,
> nanoflann, P4est, PETSC, ScalaPack, SLEPC, Sundials, Threads,
> Trilinos, Umfpack and Zlib
> and the internal version of Muparser
> Why do I get that linking error, and how can I fix it?
> Apparently it is not related to the external libraries, else
> it would show up. Is that correct?
> Thanks!
>
> -- 
> 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
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Poisson test code fails for a fe-degree >= 2

[Martin Kronbichler]

Dear Maxi,

The problem is the quadrature formula you use in 'assemble_system' and
'calculate_residual'. You should use 'QGauss
quadrature_formula(fe.degree+1)`. I am surprised that quadratic elements
even worked because the Poisson matrix for d>=2 should have some rank
deficiency with under-integration as far as I know.

Best,
Martin

On 04.04.19 11:31, 'Maxi Miller' via deal.II User Group wrote:
> I wrote a test program, based on step-5, for testing unrelated things,
> and checked the results. But when I increased the degree of the finite
> elements from 1 to 3 (or higher), the results were wrong, while the
> change from 1 to 2 just improved convergence, as expected. Did I make
> a mistake in my code here?
> -- 
> 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
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: dealii compilation on Cray machine

[Martin Kronbichler]

On 05.04.19 21:45, Bruno Turcksin wrote:

Le ven. 5 avr. 2019 à 14:48, Phani Motamarri  a écrit :

No login node is intel-broadwell and compute node is intel-skylake. Problem is 
that the supercomputer I am compiling on has no way to access compute node 
(interactive job) directly for compiling. I am not sure if dealii allows for 
cross compilation.

I've never tried that. Maybe Martin has experience with that?


Unfortunately nothing that would help - I always managed to run at least 
the configuration (cmake) on the same architecture whereas I indeed 
built (make) on a different architecture once the compilation was 
finished through and with a fixed '-march=knl' for example. But we 
should definitely support cross-compilation. I'm not sure how our 
support was - I will test this before the release.


As to the present case: Phani, have you tried to submit the compilation 
job, i.e., "cmake -D  && make -j xx" through a job file? That should 
invoke the right microarchitecture; some machines do not have all 
compile tools such as cmake available as modules on the nodes, though.


Best,
Martin

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Alternative to local_range in LinearAlgebra::distributed::Vector

[Martin Kronbichler]

> On 2/2/19 12:32 AM, Praveen C wrote:
>> Thank you for this tip. For my need, I only require the starting global 
>> index 
>> on each partition since I am using local_element to do some operations like 
>> multiplying by a local time step which is different for each cell. The 
>> local_range function gave this info. Here is a sample usage in my code
>>
>> for(; cell!=endc; ++cell)
>> if(cell->is_locally_owned())
>>        {
>> constunsignedintcell_no = cell_number (cell);
>>
>>           cell->get_dof_indices (dof_indices);
>> for(unsignedinti=0; i>           {
>> unsignedinti_loc = dof_indices[i] - local_range.first;
>>              newton_update.local_element(i_loc) = dt(cell_no) *
>>   
>> right_hand_side.local_element(i_loc) *
>>                                                   
>> inv_mass_matrix[cell_no][i];
>>           }
>>        }
>>
>> If there is a possibility, I would like to request to retain this simple 
>> function in the library.
> No :-) The reason why it was scheduled for removal is that it only works in 
> case the locally owned range is really just a single interval. But that's not 
> always the case. 
To be fair, the current implementation of LA::dist::Vector only supports
that case. However, it is only present for that particular vector type,
which is the main reason why we did not want to expose those functions.
Instead of having a list of functions exposed, we opted for only
exposing the vector's partitioner. The safest approach is to use
something like

  i_loc =
LA::dist::Vector::get_partitioner()->global_to_local(dof_indices[i]);

The partitioner is part of the interface, and global_to_local() will be
stable even if we at some point decide to support non-contiguous range
within the vector.

Best,
Martin

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Compilation with ICC18 leads to several undefined AVX/SSE-functions

[Martin Kronbichler]

Dear Maxi,

I have not yet look at all details, but what strikes my attention is 
that you added `-march=native` to your C++ compile flags. That is a gcc 
thing, so I don't know what Intel makes out of that. Can you try what 
happens if you replace it by `-xhost` which is (mostly) the Intel 
equivalent?


Another thing that I still do not understand is why this passed the 
configuration step as an AVX compile (DEAL_II_VECTORIZATION_LEVEL = 2), 
since it only seems that SSE2 is available. It might be that there is 
something left from an old installation and you need to delete 
`CMakeCache.txt` and `CMakeFiles/`.


Best,
Martin

On 01.02.19 09:29, 'Maxi Miller' via deal.II User Group wrote:
I tried to compile deal.II using ICC 18, but during compilation I got 
several errors, such as

|
~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1489):error:nooperator"+="matches 
these operands

            operand types are:__m256d +=const__m256d

~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1504):error:nooperator"-="matches 
these operands

            operand types are:__m256d -=const__m256d

~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1518):error:nooperator"*="matches 
these operands

            operand types are:__m256d *=const__m256d

~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1533):error:nooperator"/="matches 
these operands

            operand types are:__m256d /=const__m256d

~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1549):error:identifier 
"_mm256_loadu_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1562):error:identifier 
"_mm256_storeu_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1574):error:identifier 
"_mm256_stream_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1598):error:identifier 
"_mm_loadu_ps"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1600):error:identifier 
"_mm256_i32gather_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1644):error:identifier 
"_mm256_sqrt_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1659):error:identifier 
"_mm256_set1_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1661):error:identifier 
"_mm256_andnot_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1674):error:identifier 
"_mm256_max_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1687):error:identifier 
"_mm256_min_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1728):error:identifier 
"_mm256_loadu_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1732):error:identifier 
"_mm256_permute2f128_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1736):error:identifier 
"_mm256_unpacklo_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1737):error:identifier 
"_mm256_unpackhi_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1770):error:identifier 
"_mm256_permute2f128_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1774):error:identifier 
"_mm256_unpacklo_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1775):error:identifier 
"_mm256_unpackhi_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1784):error:identifier 
"_mm256_loadu_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1784):error:identifier 
"_mm256_add_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1785):error:identifier 
"_mm256_storeu_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1795):error:identifier 
"_mm256_storeu_pd"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1832):error:identifier 
"_mm256_set1_ps"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1869):error:nooperator"+="matches 
these operands

            operand types are:__m256 +=const__m256

~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1884):error:nooperator"-="matches 
these operands

            operand types are:__m256 -=const__m256

~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1898):error:nooperator"*="matches 
these operands

            operand types are:__m256 *=const__m256

~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1913):error:nooperator"/="matches 
these operands

            operand types are:__m256 /=const__m256

~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1929):error:identifier 
"_mm256_loadu_ps"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1942):error:identifier 
"_mm256_storeu_ps"isundefined


~/Downloads_ICC/dealii/include/deal.II/base/vectorization.h(1954):error:identifier 

[deal.II] Thank you!

[Martin Kronbichler]

I would like to follow up on the announcement of this release by publicly
saying again how much we appreciate the many contributions from those who
have sent it code, bug reports, fixed grammar and typos, or have helped in
any other way. Many thanks!

The ChangeLog lists at least the following people to whom thanks are due
(though I am certain that I am forgetting someone):

  Giovanni Alzetta,
  Mathias Anselmann,
  Daniel Appel,
  Alexander Blank,
  Vishal Boddu,
  Benjamin Brands,
  Pi-Yueh Chuang,
  T. Conrad Clevenger,
  Sambit Das,
  the Ginkgo developers,
  Stefano Dominici,
  Nivesh Dommaraju,
  Marc Fehling,
  Niklas Fehn,
  Isuru Fernando,
  Andreas Fink,
  Daniel Garcia-Sanchez,
  Rene Gassmöller,
  Alexander Grayver,
  Joshua Hanophy,
  Logan Harbour,
  Graham Harper,
  Daniel Jodlbauer,
  Stefan Kaessmair,
  Eldar Khattatov,
  Alexander Knieps,
  Uwe Köcher,
  Kurt Kremitzki,
  Dustin Kumor,
  Ross Kynch,
  Damien Lebrun-Grandie,
  Jonathan Matthews,
  Stefan Meggendorfer,
  Pratik Nayak,
  Lei Qiao,
  Ce Qin,
  Reza Rastak,
  Roland Richter,
  Alberto Sartori,
  Svenja Schoeder,
  Sebastian Stark,
  Antoni Vidal,
  Jiaxin Wang,
  Yuxiang Wang,
  Zhuoran Wang.


Martin,
on behalf of the whole deal.II developer team:
Daniel Arndt, Wolfgang Bangerth, Denis Davydov, Timo Heister,
Luca Heltai, Guido Kanschat, Martin Kronbichler, Matthias Maier,
Jean-Paul Pelteret, Bruno Turcksin, and David Wells

-- 
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/5a8a2837-3c10-06aa-53ec-247b6e4735d3%40gmail.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] deal.II version 9.1 released

[Martin Kronbichler]

Version 9.1.0 of deal.II, the object-oriented finite element library
awarded the
J. H. Wilkinson Prize for Numerical Software, has been released. It is
available for free under an Open Source license from the deal.II homepage at

   https://www.dealii.org/

The major changes of this release are:

  - The support for automatic and symbolic differentiation within
deal.II has
    been completely overhauled and significantly extended. The automatic
    differentiation facilities integrate a broad spectrum of functionality,
    including support in the evaluation of finite element shape functions
    (FEValues) and associated point and tensor data structures.
Functionality
    from both the ADOL-C as well as Sacado library in various modes is
    included. For symbolic differentiation, new wrappers to the
    high-performance SymEngine library allow an alternative approach
    addressing the syntax levels of expressions to compute derivatives,
rather
    than the algorithmic approach underlying automatic differentiations.

  - Full support for hp adaptivity in parallel computations. deal.II has
    had support for hp-adaptive methods since around 2005, but not yet in
    combination with MPI. In order to gain support for parallel hp
adaptivity
    a number of algorithmic issues were addressed in this release. Details
    can be found in the release paper.

  - A new HDF5 interface has been added in order to leverage
    its capabilities in terms of high-performance I/O for large amount of
    data. The HDF5 interfaces in deal.II support read and write operations
    both in serial and with MPI.

  - GPU support was significantly extended for the current release:
    Added features include preconditioners for CUDAWrappers::SparseMatrix
    objects, support for MPI-parallel CUDA data structures, and support for
    constraints in the matrix-free framework (and thus also support for
    adaptively refined meshes).

  - Four new tutorial programs have been added, step-61 demonstrating an
    implementation of the so-called weak Galerkin method, step-62
solving the
    elastic wave equation with perfectly matched layers to calculate the
    resonance frequency and bandgap of photonic crystals, step-63 on
block and
    point smoothers for geometric multigrid in the context of
    convection-diffusion problems, and step-64 implementing a Helmholtz
solver
    with matrix-free methods running on GPUs. Furthermore, a new code
gallery
    program solving a Bayesian inverse problem has been added.

  - The release contains performance improvements and bug fixes of the
    matrix-free framework and related geometric multigrid solvers. In
    particular, the implementation of the Chebyshev iteration, an often used
    smoother in the matrix-free context, has been revised to reduce vector
    accesses. Altogether, matrix-free multigrid solvers run up to 15% faster
    than in the previous version.

  - Various variants of geometric multigrid solvers and matrix-free
    implementations were run on up to 304,128 MPI ranks during the
acceptance
    phase of the SuperMUC-NG supercomputer, verifying the scalability of our
    implementations to this scale. Some geometric multigrid data structures
    were revised to avoid bottlenecks showing up with more than 100k ranks.

  - The FE_BernardiRaugel class implements the non-standard
    Bernardi-Raugel element that can be used to construct a stable
    velocity-pressure pair for the Stokes equation. The Bernardi-Raugel
    element is an enriched version of the Q_1^d element with added bubble
    functions on each edge (in 2d) or face (in 3d). It addresses the fact
    that the Q_1^d - Q_0 combination is not inf-sup stable (requiring a
    larger velocity space), and that the Q_2^d - Q_1 combination is
    stable but sub-optimal since the velocity space is too large relative to
    the pressure space to provide additional accuracy commensurate with the
    cost of the large number of velocity unknowns. The Bernardi-Raugel
    space is intermediate to the Q_1^d and Q_2^d spaces.

  - The FE_NedelecSZ class is a new implementation of the Nédélec element
    on quadrilaterals and hexahedra. It overcomes the sign conflict issues
    present in traditional Nédélec elements that arise from the edge
    and face parameterizations used in the basis functions. Therefore, this
    element should provide consistent results for general quadrilateral and
    hexahedral elements for which the relative orientations of edges and
    faces (as seen from all adjacent cells) are often difficult to
establish.

  - A new class ParsedConvergenceTable has been introduced. The class
    simplifies the construction of convergence tables, reading the options
    for the generation of the table from a parameter file. It provides a
    series of methods that can be used to compute the error given a
    reference exact solution, or the difference between two numerical
    solutions, or any other 

Re: [deal.II] Installation error on Haswell nodes on Cori at NERSC (failed AVX512 support)

[Martin Kronbichler]

Dear Steve,

>From what I can see the failure is for the expand_instantiations script
of deal.II, which is compiled as part of deal.II. It uses slightly
different flags as the full install, but assuming that you either passed
-xHASWELL or -xCORE-AVX2 to CMAKE_CXX_FLAGS it should not generate that
code. Before we look into the flags used for compilation, a basic
question: Did you start with a clean build directory without any file
left over from a more advanced architecture with AVX512 support?

Best,
Martin

On 10.07.19 15:56, Stephen DeWitt wrote:
> Hello,
> I'm trying to install deal.II on the Haswell nodes on Cori at NERSC
> using the Intel compiler. I'm using deal.II version 9.0, because
> support for a few of the function calls I make were dropped in v9.1
> and I haven't had a chance to modify the sections of the code. In my
> CMake command, I'm adding either the -xHASWELL or -xCORE-AVX2 flags to
> get the right level of vectorization for this architecture (AVX). The
> CMake output is what I expect:
>
> |
> -- Performing Test DEAL_II_HAVE_SSE2
> -- Performing Test DEAL_II_HAVE_SSE2 - Success
> -- Performing Test DEAL_II_HAVE_AVX
> -- Performing Test DEAL_II_HAVE_AVX - Success
> -- Performing Test DEAL_II_HAVE_AVX512
> -- Performing Test DEAL_II_HAVE_AVX512 - Failed
> -- Performing Test DEAL_II_HAVE_OPENMP_SIMD
> -- Performing Test DEAL_II_HAVE_OPENMP_SIMD - Success
> |
>
> However, when I try to compile I get this error:
>
> |
> [ 32%] Built target obj_umfpack_DL_STORE_release
> [ 34%] Built target obj_amd_global_release
> [ 35%] Built target obj_amd_long_release
> [ 36%] Built target obj_amd_int_release
> [ 37%] Built target obj_muparser_release
> [ 37%] Built target obj_sundials_inst
> Scanning dependencies of target obj_sundials_release
> [ 37%] Building CXX object
> source/sundials/CMakeFiles/obj_sundials_release.dir/arkode.cc.o
> [ 37%] Building CXX object
> source/sundials/CMakeFiles/obj_sundials_release.dir/ida.cc.o
> [ 37%] Building CXX object
> source/sundials/CMakeFiles/obj_sundials_release.dir/copy.cc.o
> [ 38%] Building CXX object
> source/sundials/CMakeFiles/obj_sundials_release.dir/kinsol.cc.o
> [ 38%] Built target obj_sundials_release
> [ 38%] Generating data_out_dof_data.inst
>
> Please verify that both the operating system and the processor support
> Intel(R) AVX512F, ADX, AVX512ER, AVX512PF and AVX512CD instructions.
>
> source/numerics/CMakeFiles/obj_numerics_inst.dir/build.make:91: recipe
> for target 'source/numerics/data_out_dof_data.inst' failed
> make[2]: *** [source/numerics/data_out_dof_data.inst] Error 1
> CMakeFiles/Makefile2:1860: recipe for target
> 'source/numerics/CMakeFiles/obj_numerics_inst.dir/all' failed
> make[1]: *** [source/numerics/CMakeFiles/obj_numerics_inst.dir/all]
> Error 2
> Makefile:129: recipe for target 'all' failed
> make: *** [all] Error 2
> |
>
> It seems to me that it is looking for AVX512 support and generating an
> error when it doesn't find it. But why would it look for AVX512
> support if DEAL_II_HAVE_AVX512 failed? I haven't had this issue when
> installing on other machines that support AVX but not AVX512.
>
> Thanks for any insight you have,
> Steve
> -- 
> 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/b36fd46d-2dff-41f4-8762-8ffad6ad4da1%40googlegroups.com
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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/daa0d6eb-090a-9876-bd94-c70876378752%40gmail.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] DG Matrix-free: getting face iterator form

[Martin Kronbichler]

Dear Michal,

We do not offer convenience access for all data structures on faces yet,
so feel free to suggest what you need.

That said, for the access of rho and mu on an element, we provide the
following access function:

https://www.dealii.org/current/doxygen/deal.II/classFEEvaluationBase.html#a442f1ee55e6e80a58e19d53c7184c1a7

This FEEvaluationBase::read_cell_data function expects an
AlignedVector> as the argument with as many
components as there are cell batches (plus ghosted ones for the MPI
case). So in your code, you would fill the rho and mu fields by a cell
loop (for MPI, please run it for "n_cells = this->data->n_macro_cells()
+ this->data->n_ghost_cell_batches()". Then you can access that data
both on cells and faces, i.e., with FEEvaluation and FEFaceEvaluation.

The implementation resolves the `interior_cells` and `exterior_cells`
fields for you.

Retrieving the `user_index` is a bit more involved as we do not provide
access to a face iterator.

You could do it approximately as follows (sitting on a face like in your
code below):

const unsigned int interior_cell = face2cell_info.cells_interior[element];
const Triangulation::cell_iterator my_cell =
  this->data->get_cell_iterator(interior_cell /
VectorizedArray::n_array_elements,
 interior_cell %
VectorizedArray::n_array_elements);
const Triangulation::face_iterator face_it =
my_cell->face(face2cell_info.interior_face_no);

Can you try if that works? If yes, we should add some access functions
to MatrixFree. Let us know if you are interested in helping with a pull
request.

Best,
Martin

On 05.08.19 19:34, Michał Wichrowski wrote:
> Dear all,
> I'm trying to use matrix-free with discontinuous Galerkin elements for
> variable coefficient problem.  The face integral require access to
> value of coefficient. For the cell I've done the following thing:
> //Vectors of values, stored in object:
>   AlignedVector< VectorizedArray > mu;
>   AlignedVector< VectorizedArray > rho;
> //Filling the vectors:
>   const unsigned int n_cells = this->data->n_macro_cells();
>   rho.resize  (n_cells);
>   mu.resize  (n_cells);
>   for (unsigned int cell=0; cell     {
>   for(unsigned int element=0;
> elementdata->n_components_filled(cell) ; ++element){
>   const unsigned int  ID=this->data->get_cell_iterator(cell,element,0
> )->material_id();
>   rho[cell][element] = rho_map.at(ID) ;
>   mu[cell][element]  = miu_map.at(ID) ;
>   }
>     }
> For the faces I am trying to do something like that:
> //Vectors of values, stored in object:
>   AlignedVector< VectorizedArray > mu_inner_faces;
>
> //Filling the vectors:
>   const unsigned int n_faces=this->data->n_inner_face_batches();
>   mu_inner_faces.resize(n_faces);
>   for (unsigned int face=0; face   const internal::MatrixFreeFunctions::FaceToCellTopology
>   < VectorizedArray< Number >::n_array_elements >
> face2cell_info=this->data->get_face_info(face);
>
>   for(unsigned int element=0;
>   elementdata->n_active_entries_per_face_batch(face) ;
>   ++element){
>   const unsigned int interior_cell =
> face2cell_info.cells_interior[element];
>     const unsigned int exterior_cell =
> face2cell_info.cells_exterior[element];
>   
>   ??
>
>   }
>   }
>
> If I understood correctly interior_cell and exterior_cell will be
> numbers of cell.
>  I have no idea how to get cell iterators from that information. I
> even found that MAtrixFree class contains attribute
> std::vector< std::pair< unsigned int, unsigned int > > 
> cell_level_index
> 
>
>
> but unfortunately it is private and  there is no  getter for that. 
>
> I will also need to get user_index from each face and create from that
> another AlignedVector of parameters. This operations will be done only
> once so I do not care much about efficiency.
> Best, 
> Michal
> -- 
> 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/b3c0fccd-7a76-4d90-a573-74719d281c37%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 

Re: [deal.II] DG Matrix-free: getting face iterator form

[Martin Kronbichler]

Dear Michal,

Great. I've opened issue https://github.com/dealii/dealii/issues/8460
for that to keep track of it. We won't get at it immediately but I hope
to do so in the coming weeks.

Best,
Martin

-- 
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/9073f625-871e-ddba-d642-69e38c6839c0%40gmail.com.

Re: [deal.II] Matrix-free-method for an operator depending on input from solver

[Martin Kronbichler]

Dear Maxi,

There is not really a simple answer to your request: It depends on how
good a preconditioner you need. Matrix-free methods work best if you can
use simple preconditioners in the sense that they only need some matrix
entries (if any) and the matrix-vector product (or operator evaluation
in the nonlinear case). Note that I also consider multigrid with
Chebyshev smoothing as in step-37 and step-59 a simple preconditioner.
It is up to you to find out whether your problem admits such a simple
form and iterative solvers, such that GMRES, some nonlinear
Newton-Krylov or nonlinear Krylov variant, can cope with that and not
use too many iterations. We as a mailing list cannot help with the
preconditioner because this is a research question that many people
would like an answer of. My experience is that matrix-free setups are
mostly done for performance reasons (because the mat-vec is cheaper) for
well-understood systems, like the Laplacian with multigrid or some
linear/nonlinear time-dependent equations with dominating mass matrix
terms where simple iterative methods suffice. Sometimes matrix-free is
done because one can fit larger problems in memory. But all cases where
I know it is used are done when one at least has an expectation of what
the matrix and its eigenvalue spectrum look like. If you do not have
that, because you have some difficult linear/nonlinear systems,
matrix-free setups might help provide you with a black-box view - but
they make the preconditioner selection harder, not easier, because you
restrict yourself in the choice of preconditioner to the more basic
variants.

Now to the technical part: If you need the diagonal of a matrix and rely
on an accurate representation, you will need to build the linearization
in some way or the other. This is also what step-37 or step-59 do - they
compute the local stiffness matrix, throw away everything but the
diagonal, and use that inside some solver. Tor the diagonal you would
run through all unit vectors with v, i.e., you add (1,0,0,...),
(0,1,0,...), and so on. Then you only keep the one entry where v was
one. Your approach to do this globally is certainly not going to scale
because it is an O(n^2) operation to check every unit vector. In
principle, you could do it locally element-by-element similar to step-37
also for nonlinear operators. The question is whether you have the full
action inside a single "local_apply" function such that you could start
adding unit vectors there. Given a nonlinear field u that you linearize
about, you could perturb it with e*v_i and compute the local diagonal
according to your formula. From what I understand, you have the operator
in an even more abstract way, so collecting everything in one local call
with matrix-free facilities would also be non-trivial.

Best,
Martin

On 24.07.19 15:01, 'Maxi Miller' via deal.II User Group wrote:
> I am currently trying to implement the JFNK-method in the matrix-free
> framework (by following step-37 and step-59) for solving nonlinear
> equations of the form F(u)=0. The method itself replaces the
> multiplication of the explicit jacobian J with the krylov vector v
> during solving with the approximation (F(u+e*v)-F(u))/e, which removes
> the matrix.
> The initial test using a LinearOperator worked, where I modified the
> vmult-function accordingly (c.f. my earlier questions here), but still
> I had to form the jacobian at least once (which I wanted to avoid),
> thus the test of the matrix-free framework.
> In theory I just should have to provide the vmult-function, but based
> on step-37 and step-59 I also have to provide the full operator at
> least once during initialization, as f.ex. to form the inverse
> diagonal in step-37. This would be the same as calculating
> (F(u_i+e*v_i)-F(u_i))/e, but I do not know the value of v_i.
> Thus, what would be the best approach of implementing this operator
> using the MatrixFree-framework in deal.II, without having to form a
> single matrix, and without knowing the exact value of F? Ideally it
> can be used for a preconditioner too, such as the GMG in step-37 and
> step-59.
> Thanks!
> -- 
> 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/257bb96a-ba19-48cf-a7c4-49e6f178e417%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 

Re: [deal.II] Create a LinearAlgebra::distributed::Vector from a LinearAlgebra::distributed::Vector

[Martin Kronbichler]

Dear Maxi,

Did you try to call 'operator=' in the same line as you declare the 
variable? In that case, you still go through the (missing) constructor. 
It should work if you first declare the float vector in one statement 
and then assign in the next. We have this function:


https://dealii.org/current/doxygen/deal.II/classLinearAlgebra_1_1distributed_1_1Vector.html#ad364c47c5b3b7b991490cd626c729f0b

Furthermore, there is also a function 
LinearAlgebra::distributed::Vector::copy_locally_owned_subrange_from() 
which also works with different number types, say double to float:


https://dealii.org/current/doxygen/deal.II/classLinearAlgebra_1_1distributed_1_1Vector.html#ab792ddb04b95a220e489f2d7f9eee990

(In that case, no data is exchanged, and the vectors need to have 
compatible locally owned dofs before calling.)


Best,
Martin

On 23.09.19 17:07, 'Maxi Miller' via deal.II User Group wrote:

When doing that in my code, I get the error
|
error:conversion from‘Vector’to non-scalar type 
‘Vector’requested

|

thus I assume I have to implement it myself...

Am Freitag, 13. September 2019 17:53:37 UTC+2 schrieb Wolfgang Bangerth:

On 9/13/19 8:39 AM, 'Maxi Miller' via deal.II User Group wrote:
> In my program I have to reuse old solutions (due to time
dependency),
> but would like to use them both as double values and float values
> (preconditioning does not require double values). Unfortunately
my old
> solution is only given as a
> LinearAlgebra::distributed::Vector(), but I have to use a
> LinearAlgebra::distributed::Vector. Direct conversion via
> LinearAlgebra::distributed::Vector
float_solution(old_solution)
> does not work. Are there other ways for conversion?

I think that the other classes all have conversion constructors or
conversion operator=. If this class doesn't have one, you probably
want
to add it to support what you need. It should be a relatively easy
patch
to write, and we'd of course be happy to take it, as always!

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/1e88472a-6dce-4a26-ada0-197a82adea63%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/32b70356-bde3-a7be-fd08-9880cc7539ae%40gmail.com.

Re: [deal.II] Is MatrixFreeOperators::LaplaceOperator::local_diagonal_cell() something wrong?

[Martin Kronbichler]

I should add that you can set a variable coefficient, but it is still 
scalar and thus the same to all components, so the code is correct. It 
would be different if we allowed tensor-valued coefficients (which we do 
not at the moment); then, one would really need to go through all 
dofs_per_cell entries.

-- 
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/23e11aee-02eb-41ff-b12a-020684d223d0%40googlegroups.com.

[deal.II] Re: how to implement MatrixFreeOperation::compute_diagonal in vector-valued problem?(which modified based on step-37 which used Multigrid method)

[Martin Kronbichler]

Dear m.,

As you have already found out, the code in the Stokes example should do 
what you need for computing the diagonal. Usually, the approach with 
FEEvaluation::begin_dof_values() is slightly simpler than the approach with 
the tensor you mentioned in the first email. We have computed vector-valued 
problems with matrix-free, including multigrid preconditioners which 
compute the diagonal in a similar way, so they do work.

One thing I would recommend to do is to compare the diagonal you get with 
the cell-by-cell approach with computing the global diagonal by applying 
the operator globally on all unit vectors, i.e., outside the matrix-free 
loop. This will only work up to a few 10s of thousands of vectors because 
it is an n^2 operation and obviously the wrong choice for production runs, 
but it is a useful debugging approach in case you do not know for sure 
whether the interior things are set up correctly (well, constraints might 
still possibly mess this up, but that is a different question).

Regarding your last question:

> is the Multigrid still work well in vector-valued problem? or in 
> multi-physics couple problem? 
> I want to use it to solve the couple equations(Allen-Cahn and mechanics 
> equalibrium equations
>

This depends a lot on the equation. For coupled problems with non-trivial 
coupling, the answer is generally "no, multigrid alone will not work out of 
the box". But that is not specific to the implementation (e.g. whether 
matrix-free or not) but due to the mathematics of the underlying equations 
and multigrid convergence theory. In general, you can expect that simple 
smoothers like Chebyshev will not like multiphysics component coupling, and 
you need to address smoothing "by hand". What people do is to iterate 
between the field in an appropriate way and use simple smoothers on each 
component of the multiphysics problem, as e.g. described in the paper here: 
https://www.sciencedirect.com/science/article/pii/S0045782516307575 . But I 
would say that the field is still pretty ad hoc. Then, the next question 
with Allan-Cahn is whether something like Chebyshev+point-Jacobi still 
works well in the nonlinear case with the typical contrasts in the terms. 
Maybe yes, maybe no. The consensus in the field is that for strongly 
varying coefficients you need to either construct very specific smoothers, 
or, more commonly, use operator-dependent level transfer operators that are 
able to react to those case. AMG is supposed to handle some cases better 
than GMG with naive transfer as we provide in deal.II - but AMG is far from 
a black box. I have done work with Cahn-Hilliard, which involves a block 
system of two components, and there you definitely do not want to simply 
apply multigrid to the 2x2 block system as a whole, but use some algebraic 
manipulations on the blocks with a basic preconditioner of some related 
operator instead: 
https://www.sciencedirect.com/science/article/pii/S0898122112004191 .

I hope this helps.

Best,
Martin

-- 
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/04cca9d8-cda3-4f0a-9ffd-7c67cf8063b3%40googlegroups.com.

Re: [deal.II] Is MatrixFreeOperators::LaplaceOperator::local_diagonal_cell() something wrong?

[Martin Kronbichler]

> In the first loop inside the cell loop index i go through from 0 to
> phi.dofs_per_component (highlighted) as well as the j index,
> but the function phi.begin_dof_values()return a pointer whose index go
> throught from 0 the phi.dofs_per_cell, that is said
> *in n_component > 1 case the index after phi.dofs_per_component never
> been touched( the total is dofs_per_compoenent*n_component).*
> I have found an old discussion
> in 
> https://groups.google.com/forum/#!searchin/dealii/MatrixFreeOperators$3A$3ALaplaceOperator%7Csort:date/dealii/k133wjcKM70/Jw-reYrsAwAJ\
> but things change now.

The code here should be correct, because we are looking at a
constant-coefficient Laplacian here, where all entries to the diagonal
are equal for the various components. So even though the local
computation does indeed compute all contributions (but with some
spurious dummy entries in the components higher than zero; which
eventually are never used), but only get the zero-th component of it and
duplicate it over all entries. This is done to save some work, given
that we know that the entry is the same; in case you do an operator that
couples between the components, one of course has to go over all
phi.dofs_per_cell components both when zeroing entries and in the loop
over the columns i.

Best,
Martin

-- 
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/98b32fdc-2208-0745-fbc6-5420975ed08f%40gmail.com.

Re: [deal.II] Extracting DoF indice from FEEvaluation / Distributing from FEEvaluation to row of a matrix

[Martin Kronbichler]

Dear Michal,

To answer the first request, we can provide an iterator from 
`MatrixFree` via `MatrixFree::get_cell_iterator(cell_index, v)`:


https://www.dealii.org/developer/doxygen/deal.II/classMatrixFree.html#adb12e97b3905c1cb12d6f90527604476

That would work for the cell part. However, we are currently lacking an 
equivalent functionality for the faces (it should be easy to add, though).


Regarding the second request to let `FEEvaluation` handle the transfer 
into a matrix: Peter Munch already did some work into that direction. 
I'm adding him in CC so we can discuss together how to best provide a 
solution for your problem. In the meantime, I'm adding an issue as well 
to get that information.


Best,
Martin

On 03.10.19 21:01, Michał Wichrowski wrote:

Dear all,
I need to construct a sparse matrix from MatrixFree operator. I want 
to do it in a similar way to present in step-37 for computing 
diagonal. The only difference will be that the computed values from  
FEEValuation and FEFaceEvaluation have to be distributed to the 
specified row of a matrix. The code will look like, I need to get 
local_dof_indices somehow:


// Loop is executed without threads
for (unsigned int cell=0; cellphi.submit_dof_value(VectorizedArray 
(), 
j);

phi.submit_dof_value(make_vectorized_array(1.), i);
// local integration operator:
integrate_cell(phi,cell);
phi.integrate (false, true);
// Here I have phi filled with a row of the matrix
// now I need to know global indices of dofs...
// And I have no clue how to obtain local_dof_indices needed below.
for (unsigned int j=0; jdata->n_components_filled(cell); ++v)
matrix(local_dof_indices[i][v],local_dof_indices[j][v] 
)+=phi.get_dof_value(j)[v];

}
}
}
The code for FEFaceVAlues will be very similar.
Moreover, I'm working with a vector-valued problem so 
local_dof_indices will need 3 indices ( local dof number, index in 
vectorization and component of the vector)
It would be great it FEEvaluation had a function for distributing to a 
matrix so that two loops are replaced with 
phi.distribute_local_to_global (matrix, i);
I only need to do that somehow, it does not have to be very efficient. 
There are three main purposes for having that function:

1) coarse matrix for direct solver
2) initialization of preconditioner based on sparse matrix
3) system matrix for validating code
Having this functionality would allow me to write a class for an 
arbitrary matrix-free operator that form local integration formula 
will be able to automatically:

1) compute matrix-vector product
2) compute diagonal
3) compute sparse matrix of operator.
The operator will be specyfied by defining a function :
void integrate_cell(FEEValuation & phi, unsigned int cell);
And similar ones for face and boundary in case of DG. I think it will 
be possible for both vector and scalar based 2nd order operators. Also 
non-square operator could be possible. And using this, the block 
operators may be constructed (eg. Stokes).

Best,
Michał
--
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/e3f68a4b-4361-4167-9514-d89e746bfd8b%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/a957cc63-bf0a-b16e-95db-1373f60ee537%40gmail.com.

Re: [deal.II] Matrix-free DG with adaptive refinements

[Martin Kronbichler]

Dear Michal,

Thanks for reporting this problem - this is indeed a bug. It seems we 
did not think carefully about the case with periodic boundary conditions 
and hanging nodes. I have opened an issue at deal.II so we can try to 
fix it soon: https://github.com/dealii/dealii/issues/8905


Best,
Martin

On 12.10.19 17:29, Michał Wichrowski wrote:

Dear all,
I modified step-59 function run() so that in only makes grid and does 
local refinement and initializes MatrixFree object.  This triggers an 
assertion when calling

        mg_mf_storage_level->reinit(dof_handler,
                                    dummy,
QGauss<1>(fe.degree + 1),
                                    additional_data);
in system_setup().

Function run():
            Point upper_right;
            upper_right[0] = 2.5;
            for (unsigned int d = 1; d < dim; ++d)
              upper_right[d] = 2.8;
            GridGenerator::hyper_rectangle(triangulation,
 Point(),
                                           upper_right);
triangulation.begin_active()->face(0)->set_boundary_id(10);
triangulation.begin_active()->face(1)->set_boundary_id(11);
triangulation.begin_active()->face(2)->set_boundary_id(0);
            for (unsigned int f = 3; f < 
GeometryInfo::faces_per_cell; ++f)

triangulation.begin_active()->face(f)->set_boundary_id(1);

std::vector::cell_iterator>>
              periodic_faces;
            GridTools::collect_periodic_faces(
              triangulation, 10, 11, 0, periodic_faces);
triangulation.add_periodicity(periodic_faces);

            triangulation.refine_global(1);

            Vector 
estimated_error_per_cell(triangulation.n_active_cells());
            for (unsigned int i = 0; i < 
estimated_error_per_cell.size(); ++i)

              estimated_error_per_cell[i] = std::sin(i) * std::sin(i);

parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
                triangulation, estimated_error_per_cell, 0.3, 0.0);
triangulation.execute_coarsening_and_refinement();
        setup_system();

Error message:


An error occurred in line <4993> of file 
 in function
    void dealii::FESubfaceValues::reinit(const typename 
dealii::Triangulation::cell_iterator&, unsigned int, 
unsigned int) [with int dim = 3; int spacedim = 3; typename 
dealii::Triangulation::cell_iterator = 
dealii::TriaIterator >]

The violated condition was:
    subface_no < cell->face(face_no)->n_children()
Additional information:
    Index 0 is not in the half-open range [0,0). In the current case, 
this half-open range is in fact empty, suggesting that you are 
accessing an element of an empty collection such as a vector that has 
not been set to the correct size.


Stacktrace:
---

I have no idea what may be causing this.
--
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/69586e52-7e99-41e4-a7f1-930cffd6c046%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/556cf82f-79bf-6b3b-3265-4e84501f78f1%40gmail.com.

Re: [deal.II] Facilities for computing the "goodness" of a mesh (i.e. something like a mesh-quality indicator)

[Martin Kronbichler]

Dear Krishna,

We have recently added one such function, namely 
GridTools::compute_aspect_ratio_of_cells(),


https://dealii.org/developer/doxygen/deal.II/namespaceGridTools.html#a2a9d8801d952ea8afa9ec1573229a329

Best,
Martin

On 15.01.20 18:26, Krishnakumar Gopalakrishnan wrote:


In the results/extension section of Step-6 as well in Step-49, there 
is a great deal of discussion about obtaining good meshes are for a 
problem, with the caveat that this is a highly specialized topic. 
Step-49 is full of warnings about things that can go wrong when 
applying certain transformations, and it seems to implant doubts in my 
mind that my meshes are never going to be correct :)



Given this, in the spirit of the Kelly Indicator for solution errors 
in a cell, is there a similar single-metric for the quality/goodness 
of each cell in a mesh? Maybe there is something based on the smallest 
internal angle (to prevent slivers, for Jacobian of the mapping to be 
invertible, to prevent inverted cells etc).



After all, error estimation is a complicated topic. nevertheless the 
Kelly Indicator seems to do a reasonable job for quite a few 
equations. If currently not implemented in Deal.II, is there any 
journal article that discusses research in this direction?




Krishna

--
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/3602a53a-a98a-4a5c-a512-a6d179c8dcec%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/36574d22-6fa9-e153-637a-2641e798756d%40gmail.com.

Re: [deal.II] Segfault when calling gather_evaluate() or read_dof_values_plain()

[Martin Kronbichler]

Dear Maxi,

It looks like you are using a scalar finite element (FE_Q), but you 
set FEEvaluation to operate on a vectorial field, i.e.,


        FEEvaluation phi(data);

Notice the 4th template parameter "dim", which should be one. I agree it 
is unfortunate that we do not provide a better error message, so I 
opened an issue for it, https://github.com/dealii/dealii/issues/9312


If you switch the parameter to 1, it should work. However, I want to 
point out that you use a continuous element with 
CellwiseInverseMassMatrix - this will not give you a valid matrix 
inverse. You need to either use a diagonal mass matrix (see step-48) or 
use the cellwise inverse as a preconditioner (when combined with 
suitable scaling) and solve the mass matrix iteratively.


Best,
Martin

On 14.01.20 10:49, 'Maxi Miller' via deal.II User Group wrote:

I could narrow it down to the function
|
// same as above for has_partitioners_are_compatible == true
template<
typenameVectorType,
typenamestd::enable_if::value,
VectorType>::type *=nullptr>
inlinevoid
  check_vector_compatibility(
constVectorType&                     vec,
constinternal::MatrixFreeFunctions::DoFInfo_info)
{
(void)vec;
(void)dof_info;
Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
ExcMessage(
"The parallel layout of the given vector is not "
"compatible with the parallel partitioning in MatrixFree. "
"Use MatrixFree::initialize_dof_vector to get a "
"compatible vector."));
}
|

located in vector_access_internal.h. Apparently I have a segfault in 
the function

|
template
inlinebool
Vector::partitioners_are_compatible(
constUtilities::MPI::Partitioner)const
{
returnpartitioner->is_compatible(part);
}
|

located in la_parallel_vector.templates.h, and thereby directly 
stopping the program without triggering the assert()-function. The 
direct gdb output for this function is

|
#3  0x700c0c5c in 
dealii::LinearAlgebra::distributed::Vectordealii::MemorySpace::Host>::partitioners_are_compatible (this=0x0, 
part=...) at 
~Downloads/git-files/dealii/include/deal.II/lac/la_parallel_vector.templates.h:2021

|

I am not sure what that means. Could it point to a nullpointer 
somewhere (after this=0x0)? Though, when putting a breakpoint to the 
first function and printing the involved variables there, I get

|
(gdb)printdof_info.vector_partitioner
$2 =std::shared_ptr(usecount 
3,weak count 0)={get()=0x7a60c0}

(gdb)printvec.partitioner
$3 =std::shared_ptr(usecount 
3,weak count 0)={get()=0x7a60c0}

|

i.e. no null ptr. Could the problem be there, or somewhere else?
Thanks!

Am Montag, 13. Januar 2020 21:41:54 UTC+1 schrieb Wolfgang Bangerth:

On 1/13/20 8:41 AM, 'Maxi Miller' via deal.II User Group wrote:
> Therefore, I assume that I have some bugs in my code, but am not
experienced
> enough in writing matrix-free code for being able to debug it
myself. What
> kind of approach could I do now, or is there simply a glaring
bug which I did
> not see yet?

I haven't even looked at the code yet, but for segfaults there is
a relatively
simple cure: Run your code in a debugger. A segmentation fault is
a problem
where some piece of code tries to access data at an address
(typically through
a pointer) that is not readable or writable. The general solution
to finding
out what the issue is is to run the code in a debugger -- the
debugger will
then stop at the place where the problem happens, and you can
inspect all of
the values of the pointers and variables that are live at that
location. Once
you see what variable is at fault, the problem is either already
obvious, or
you can start to think about how a pointer got the value it has
and why.

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/d55e53b2-a0c8-4cc0-ab84-7521f8b23d81%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 

Re: [deal.II] MtrixFree for biharmonic problem: submit_laplacian

[Martin Kronbichler]

Dear Michal,

That is great news. There is interest in this functionality also by 
Julius Witte (I added him here) and there were some steps done.


The main difficulty in terms of user interface is the fact that the 
submit_laplacian() is a bit nastier than the other functions: on 
non-affine geometries, the Laplacian (or Hessian) of the test function 
in real space gives contributions to both the Laplacian and the gradient 
in reference space. So we either need a `submit_laplacian()` function 
that both inserts the desired contribution to test by the gradient and 
the Laplacian - which breaks symmetry with the 
get_gradient()/get_laplacian() methods - or we need to have a second 
field that stores the accumulated part for the reference-cell gradient 
of the test function and merges this inside 'integrate'.


When you say "we already have most functionality needed to to it (DG, 
computing laplacian on both faces and cells)", do you imply that you 
have already started with the work needed for get_laplacian() on faces? 
At least on dealii master, we lack both the ability to evaluate the 
second derivatives in tangential directions in reference space and the 
Jacobian gradient needed for non-affine cell shapes.


In any case, I would be happy to help and join the efforts, and I think 
so is Julius. We have an imminent release in a few weeks for deal.II 
(scheduled feature freeze May 8), so the goal would be to try to achieve 
this soon after the release.


Best,
Martin

On 25.04.20 12:28, Michał Wichrowski wrote:

Dear all,

I was trying to implement MatrixFree biharmonic solver. We already 
have most of functionality needed to do it (DG, computing laplacian on 
both faces and cells), but the one thing missing is extending 
FEEvaluation and FEFaceEvaluation to support integration of 2nd order 
formulae: submit_laplacian() and extending integrate() would resolve 
the problem and allow matrix-free version  of step-47.


The biharmonic solver may be also useful in other parts of library, 
especially handling mesh deformation and transforms.


Best,
Michał
--
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/2448bf19-ee2f-4bdb-8ddd-f192d064db0e%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/2657163a-3aa2-5481-d3e1-09775ba17706%40gmail.com.

Re: [deal.II] extract_dofs() extremely slow in for single-thread run?

[Martin Kronbichler]

Indeed the other function was fast. I looked into the function already 
and have a fix posted:


https://github.com/dealii/dealii/pull/10883

Best,
Martin

--
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/377979f8-d10e-065e-f157-abecadebbe92%40gmail.com.

Re: [deal.II] Geometric Conservation Law

[Martin Kronbichler]

Dear Alex,

This has been on my list of things to implement and verify with deal.II 
over a range of examples for quite a while, so I'm glad you bringing the 
topic up. It is definitely true that our way to define Jacobians does 
not take those identities into account, but I believe we should add 
support for them. The nice thing is that only some local computations 
are necessary, so having the option to use it in the polynomial mapping 
classes would be great. If you would be interested in this feature and 
trying to implement things, I'd be happy to guide you to the right 
places in the code.


Best,
Martin

On 17.06.20 06:02, Alexander Cicchino wrote:

Thank you for responding Wolfgang Bangerth.

The GCL condition comes from the discretized scheme satisfying 
free-stream preservation. I will demonstrate this for 2D below, (can 
be interpreted for spectral, DG, finite difference, finite volume etc):
Consider the conservation law: \frac{\partial W}{\partial t} + 
\frac{\partial F}{\partial x} +\frac{\partial G}{\partial y} =0

Transforming this to the reference computational space (x,y)->(\xi, \eta):
J*\frac{\partial W}{\partial t} + J*\frac{ \partial \xi}{\partial x} * 
\frac{\partial F}{\partial \xi} + J * \frac{ \partial \eta}{\partial 
x}* \frac{\partial F}{\partial \eta} + J * \frac{ \partial 
\xi}{\partial y} * \frac{\partial G}{\partial \xi} + J*\frac{ \partial 
\eta}{\partial y}*\frac{\partial G}{\partial \eta}

Putting this in conservative form results in:
J\frac{\partial W}{\partial t} + \frac{\partial}{\partial \xi} ( 
J*F*\frac{\partial \xi}{\partial x} +J*G*\frac{\partial \xi}{\partial 
y} ) + \frac{\partial}{\partial \eta} ( J*F*\frac{\partial 
\eta}{\partial x} +J*G*\frac{\partial \eta}{\partial y} ) - F*( GCL in 
x) - G*(GCL in y) =0


where GCL in x = \frac{\partial }{\partial \xi} ( det(J)* 
\frac{\partial \xi

 }{\partial x}) + \frac{\partial }{\partial \eta}( det(J)* \frac{\partial
\eta}{\partial x} )
similarly for y.

So for the conservative numerical scheme to satisfy free stream 
preservation, the GCL conditions must go to zero.
For linear grids, there are no issues with the classical definition 
for the inverse of the Jacobian, but what Kopriva had shown (before 
him Thomas and Lombard), was that the metric Jacobian has to be 
calculated in either a "conservative curl form" or an "invariant curl 
form" since it reduces the GCL condition to the divergence of a curl, 
which is always discretely satisfied. In the paper by Kopriva, he 
shows this, an example in 3D:

 Analytically
J*\frac{\partial \xi}{\partial x} = \frac{\partial z}{\partial \zeta} 
* \frac{\partial y}{\partial \eta} - \frac{\partial z}{\partial \eta} 
* \frac{\partial y}{\partial \zeta}


but the primer doesn't satisfy free-stream preservation while the 
latter ("conservative curl form") does.


I will put together a unit test for a curvilinear grid.

Thank you,
Alex

On Tuesday, June 16, 2020 at 10:24:59 PM UTC-4, Wolfgang Bangerth wrote:


Alexander,

> I am wondering if anybody has also found that the inverse of the
Jacobian from
> FE Values, with MappingQGeneric does not satisfy the Geometric
Conservation
> Law (GCL), in the sense of:
>
> Kopriva, David A. "Metric identities and the discontinuous
spectral element
> method on curvilinear meshes." /Journal of Scientific Computing/
26.3 (2006): 301.
>
> on curvilinear elements/manifolds in 3D.
> That is:
> \frac{\partial }{\partial \hat{x}_1} *det(J)* \frac{\partial
\hat{x}_1
> }{\partial x_1} + \frac{\partial }{\partial \hat{x}_2} *det(J)*
\frac{\partial
> \hat{x}_2}{\partial x} + \frac{\partial }{\partial \hat{x}_3} *
> det(J)*\frac{\partial \hat{x}_3 }{\partial x_1} != 0 (GCL says
it should =0,
> similarly for x_2 and x_3)
>
> If so or if not, also, has anybody found a remedy to have the
inverse of the
> Jacobian from FE Values with MappingQGeneric to satisfy the GCL.

I'm not sure any of us have ever thought about it. (I haven't --
but I really
shouldn't speak for anyone else.) Can you explain what this equality
represents? Why should it hold?

I'm also unsure whether we've ever checked whether it holds
(exactly or
approximately). Can you create a small test program that
illustrates the
behavior you are seeing?

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 

Re: [deal.II] Geometric Conservation Law

[Martin Kronbichler]

Dear Alex,

Great! The first thing we need to know is the equation. I had a quick 
look in the paper by Kopriva and I think we want to use either equation 
(36) or (37), depending on whether we consider the conservative or 
invariant curl form, respectively. In either case, it appears that we 
need to do this in a two-step procedure. The first step is to compute 
X_l and \nabla_\xi X_m, which in deal.II speak are the "q_points" and 
"Jacobians". The implementation in mapping_q_generic.cc is a bit 
involved because we have a slow algorithm (working for arbitrary 
quadrature rules) and a fast one for tensor product quadrature rules. 
Let us consider the fast one because I think we have most ingredients 
available, whereas we would need to fill additional fields for the slow 
algorithm. The source code for those parts is here:


https://github.com/dealii/dealii/blob/9e05a87db802ecd073bf7567d77f3491170d84b4/source/fe/mapping_q_generic.cc#L1463-L1592

I skipped the part on the Hessians (second derivative of transformation) 
because we won't need them. The important parts here are the extractions 
of the positions in line 1554 and the extraction of the gradients 
(contravariant transformations) in line 1575. Those two parts will be 
the starting point for the second phase we need to do in addition: 
According to the algorithm by Kopriva, we need to define this as the 
curl of the discrete interpolation of X_l \nabla_\xi X_m. To get the 
curl, we need another round through the SelectEvaluator::evaluate() call 
in that function to get the reference-cell gradient of that object, from 
which we can then collect the entries of the curl. To call into evaluate 
one more time, we also need a new data.shape_info object that does the 
collocation evaluation of derivatives. That should only be two lines 
that I can show you how and where to add, so let us not worry about that 
part now. What is important to understand first (in terms of index 
notation vs tensor notation) is the dimension of the object. I believe 
that X_l \nabla_\xi X_m is a rank-two tensor, so it has dim*dim 
components, and we compute the gradient that gives us a dim * dim * dim 
tensor. Taking the curl in the derivative and inner tensor dimension 
space, we get rid of one component and up with a dim * dim tensor as 
expected. The final step we need to do is to divide by the determinant 
of the Jacobian (contravariant vectors), because the inverse Jacobian in 
deal.II does not yet pre-multiply with the determinant.


Does that procedure sound reasonable to you? If yes, we could start 
putting together the ingredients. It would be good to have a mesh (the 
curvilinear case you were mentioning) where we can test those formulas.


Best,
Martin

On 17.06.20 18:37, Alexander Cicchino wrote:

Dear Martin,

Thank you for your response. Yes I agree that only some local 
computations are necessary to implement the identities.
Yes I would be interested in this feature and trying to implement it. 
Do you have any suggestions on where I should start and overall 
practices I should follow?


Thank you,
Alex

On Wednesday, June 17, 2020 at 1:19:29 AM UTC-4, Martin Kronbichler 
wrote:


Dear Alex,

This has been on my list of things to implement and verify with
deal.II over a range of examples for quite a while, so I'm glad
you bringing the topic up. It is definitely true that our way to
define Jacobians does not take those identities into account, but
I believe we should add support for them. The nice thing is that
only some local computations are necessary, so having the option
to use it in the polynomial mapping classes would be great. If you
would be interested in this feature and trying to implement
things, I'd be happy to guide you to the right places in the code.

Best,
Martin

On 17.06.20 06:02, Alexander Cicchino wrote:

Thank you for responding Wolfgang Bangerth.

The GCL condition comes from the discretized scheme satisfying
free-stream preservation. I will demonstrate this for 2D below,
(can be interpreted for spectral, DG, finite difference, finite
volume etc):
Consider the conservation law: \frac{\partial W}{\partial t} +
\frac{\partial F}{\partial x} +\frac{\partial G}{\partial y} =0
Transforming this to the reference computational space
(x,y)->(\xi, \eta):
J*\frac{\partial W}{\partial t} + J*\frac{ \partial \xi}{\partial
x} * \frac{\partial F}{\partial \xi} + J * \frac{ \partial
\eta}{\partial x}* \frac{\partial F}{\partial \eta} + J * \frac{
\partial \xi}{\partial y} * \frac{\partial G}{\partial \xi} +
J*\frac{ \partial \eta}{\partial y}*\frac{\partial G}{\partial \eta}
Putting this in conservative form results in:
J\frac{\partial W}{\partial t} + \frac{\partial}{\partial \xi} (
J*F*\frac{\partial \xi}{\partial x} +J*G*\frac{\partial
\xi}{\partial y} ) + \frac{\partial}{\partial \eta} (
J*F*\frac{\p

Re: [deal.II] Geometric Conservation Law

[Martin Kronbichler]

Dear Alex,

Great! I would suggest to start by simply adding new code to the 
maybe_update_q_points_Jacobians_... function with the option to turn it 
off or on. Depending on how the final implementation will look like we 
might want to move that to a separate place, but I think it will be less 
repetitive if we use a single place.


Best,
Martin

On 22.06.20 19:59, Alexander Cicchino wrote:

Dear Martin,

Thank you very much! I have been working on making the test case not 
depend on our in house flowsolver's functions.
I think that implementing Eq. 36 the "conservative curl" form would be 
sufficient.
Yes this procedure sounds perfect to me, and I agree with the 
dimension of the object described. I have been going through the 
source code that you sent to familiarize myself with the objects. 
Should I be adding to the function 
maybe_update_q_points_Jacobians_and_grads_tensor or should I create a 
new function for it?


Thank you,
Alex

On Friday, June 19, 2020 at 5:09:14 AM UTC-4, Martin Kronbichler wrote:

Dear Alex,

Great! The first thing we need to know is the equation. I had a
quick look in the paper by Kopriva and I think we want to use
either equation (36) or (37), depending on whether we consider the
conservative or invariant curl form, respectively. In either case,
it appears that we need to do this in a two-step procedure. The
first step is to compute X_l and \nabla_\xi X_m, which in deal.II
speak are the "q_points" and "Jacobians". The implementation in
mapping_q_generic.cc is a bit involved because we have a slow
algorithm (working for arbitrary quadrature rules) and a fast one
for tensor product quadrature rules. Let us consider the fast one
because I think we have most ingredients available, whereas we
would need to fill additional fields for the slow algorithm. The
source code for those parts is here:


https://github.com/dealii/dealii/blob/9e05a87db802ecd073bf7567d77f3491170d84b4/source/fe/mapping_q_generic.cc#L1463-L1592

<https://github.com/dealii/dealii/blob/9e05a87db802ecd073bf7567d77f3491170d84b4/source/fe/mapping_q_generic.cc#L1463-L1592>

I skipped the part on the Hessians (second derivative of
transformation) because we won't need them. The important parts
here are the extractions of the positions in line 1554 and the
extraction of the gradients (contravariant transformations) in
line 1575. Those two parts will be the starting point for the
second phase we need to do in addition: According to the algorithm
by Kopriva, we need to define this as the curl of the discrete
interpolation of X_l \nabla_\xi X_m. To get the curl, we need
another round through the SelectEvaluator::evaluate() call in that
function to get the reference-cell gradient of that object, from
which we can then collect the entries of the curl. To call into
evaluate one more time, we also need a new data.shape_info object
that does the collocation evaluation of derivatives. That should
only be two lines that I can show you how and where to add, so let
us not worry about that part now. What is important to understand
first (in terms of index notation vs tensor notation) is the
dimension of the object. I believe that X_l \nabla_\xi X_m is a
rank-two tensor, so it has dim*dim components, and we compute the
gradient that gives us a dim * dim * dim tensor. Taking the curl
in the derivative and inner tensor dimension space, we get rid of
one component and up with a dim * dim tensor as expected. The
final step we need to do is to divide by the determinant of the
Jacobian (contravariant vectors), because the inverse Jacobian in
deal.II does not yet pre-multiply with the determinant.

Does that procedure sound reasonable to you? If yes, we could
start putting together the ingredients. It would be good to have a
mesh (the curvilinear case you were mentioning) where we can test
those formulas.

Best,
Martin

On 17.06.20 18:37, Alexander Cicchino wrote:

Dear Martin,

Thank you for your response. Yes I agree that only some local
computations are necessary to implement the identities.
Yes I would be interested in this feature and trying to implement
it. Do you have any suggestions on where I should start and
overall practices I should follow?

Thank you,
Alex

On Wednesday, June 17, 2020 at 1:19:29 AM UTC-4, Martin
Kronbichler wrote:

Dear Alex,

This has been on my list of things to implement and verify
with deal.II over a range of examples for quite a while, so
I'm glad you bringing the topic up. It is definitely true
that our way to define Jacobians does not take those
identities into account, but I believe we should add support
for them. The nice thing is that only some local computations
are ne

Re: [deal.II] Suggestions for "cheap" distributed preconditioner for transient systems?

[Martin Kronbichler]

Dear Bruno,




Yes, I am applying ILU to the whole system, without anything else 
really. I think this is relatively traditional in GLS and VMS 
approaches (although I would not be able to say if that is a good or a 
bad thing).


If I understand correctly,  you mean formulating the problem using a 
BlockMatrix and a BlockVector, then preconditioning each matrix of the 
block independently?
I guess this would be best achieved using the linear operators within 
deal.II right?


Yes, you would apply block vectors and block matrices. Concerning the 
linear operators that would be one option, but it is not really a 
necessity. For example you can have a look at step-31 and the class 
BlockSchurPreconditioner in particular, see its implementation here: 
https://github.com/dealii/dealii/blob/2a6983945af78ee20b10757ab010040b9f2a2f89/examples/step-31/step-31.cc#L407-L417



What would be the added benefit of that? I would presume that since 
each block have more "coherent units", this would at least to a better 
conditioning of the final system, but is there anything else there?

Is there any literature or elements on this?


The benefit is that you can tailor the components to the problem at 
hand. If you have a Navier-Stokes problem with small time steps, at some 
point you will need a Poisson-type operator in the pressure variable to 
reflect the incompressibility constraint. To resolve that global 
coupling, you need some form of multigrid (or other strong 
preconditioners). However, the effect is really only needed for the 
pressure, so it might be inefficient to also apply multigrid on the 
velocity components, and instead get around much more cheaply by doing 
an ILU there or even just the diagonal. For small time steps, the 
diagonal should be a good approximation, even though it might not be the 
most efficient choice.


Regarding literature: I would recommend the book by Elman, Silvester and 
Wathen, https://doi.org/10.1093/acprof:oso/9780199678792.001.0001 , or 
the review paper by Benzi, Golub and Liesen, 
https://dx.doi.org/10.1017/S0962492904000212


In my own experience, and ILU on the velocity block and a 
Cahouet-Chabard approach for the pressure Schur complement works well 
for medium scale simulations, i.e., polynomial degrees 2 and 3 and up to 
a few billion of unknowns and maybe 10k processor cores, see also here: 
https://doi.org/10.1177/1094342016671790 . For very large sizes, ILU is 
not really that good any more as Timo said, so one can reduce cost by 
just using the mass matrix. And geometric multigrid works better on that 
scale than algebraic multigrid used in that paper. But if you want 
efficiency for large scales and with higher Reynolds numbers, it is hard 
to beat splitting-type solvers.


Best,
Martin

--
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/5d14a21a-657e-b0e2-7109-b0da017ad45a%40gmail.com.

Re: [deal.II] Pressure-correction scheme: Behaviour of the H1 norm on the Taylor-Green vortex

[Martin Kronbichler]

Dear Jose,

While I have not experimented in detail with the step-35 program, we 
have done extensive studies on similar problems in 
https://doi.org/10.1016/j.jcp.2017.09.031 (or 
https://arxiv.org/abs/1706.09252 for an earlier preprint version of the 
same manuscript) including the pressure correction scheme. While the 
spatial discretization is DG where some of the issues are , the 
experience from our experiments suggests that the pressure correction 
scheem should not behave too differently from other time discretization 
schemes with similar ingredients (say BDF-2 on the fully coupled 
scheme). In other words, I would expect second order convergence in the 
pressure for Taylor-Hood elements. I should note that there are some 
subtle issues with boundary conditions in projection schemes, so I 
cannot exclude some hidden problems with the step-35 implementation or 
the way you set up the experiments.


Best,
Martin

On 18.10.20 18:21, jose.a...@gmail.com wrote:


Hello everyone,

I am working on a parallel solver based upon the incremental 
rotational pressure-scheme  as seen in the overview paper 
 
of Guermond and implemented instep-35 
.  To 
verify the code, the following problems were replicated


  * Numerical test proposed in section 3.7.2 of Guermond's paper. It
is based on the method of manufactured solutions. The domain is
the unit square and a circle of radius 0.5. Dirichlet boundary
conditions are imposed on the whole boundary for the velocity. The
pressure field is constrained by setting its mean value to zero.

  * The Taylor-Green vortex. The domain is the square (0,1)^2. The
velocity and pressure field are constrained by periodic boundary
conditions. Furthermore, the pressure field is constrained by
setting its mean value to zero.

In order to compare the numerical pressure solution to its analytical 
counterpart, the mean value of the numerical pressure field is made to 
match that of the analytical one at the end of each time step.


The results for the Guermond's numerical tests for the square domain

                               Velocity convergence table
==
level    dt     cells     dofs       hmax           L2             
  H1 Linfty
    7 1.00e-01  8572 132098 1.10e-02 8.557567e-03     - 
6.456461e-02  -   1.103922e-02     -
    7 2.50e-02  8572 132098 1.10e-02 5.354964e-04 -2.00 5.090698e-03 
-1.83 6.822255e-04 -2.01
    7 6.25e-03  8572 132098 1.10e-02 3.358435e-05 -2.00 4.221552e-04 
-1.80 4.245914e-05 -2.00
    7 1.56e-03  8572 132098 1.10e-02 2.105764e-06 -2.00 3.632448e-05 
-1.77 2.656204e-06 -2.00


                               Pressure convergence table
==
level    dt     cells   dofs       hmax           L2           
  H1 Linfty
    7 1.00e-01  8572 16641 1.10e-02 5.543160e-03     - 2.687030e-02   
  -    3.649990e-02     -
    7 2.50e-02  8572 16641 1.10e-02 3.394874e-04 -2.01 3.684615e-03 
-1.43 4.209359e-03 -1.56
    7 6.25e-03  8572 16641 1.10e-02 2.126561e-05 -2.00 2.738788e-03 
-0.21 3.958393e-04 -1.71
    7 1.56e-03  8572 16641 1.10e-02 3.049911e-06 -1.40 2.727803e-03 
-0.00 2.751900e-05 -1.92


and for the circle domain

                               Velocity convergence table
==
level    dt     cells  dofs      hmax L2               
   H1 Linfty
    6 1.00e-01 10979 164354 1.34e-02 7.593191e-03     - 5.404824e-02   
  -    9.465067e-03     -
    6 2.50e-02 10979 164354 1.34e-02 4.879026e-04 -1.98 4.298417e-03 
-1.83 6.072627e-04 -1.98
    6 6.25e-03 10979 164354 1.34e-02 3.069737e-05 -2.00 3.539538e-04 
-1.80 3.808221e-05 -2.00
    6 1.56e-03 10979 164354 1.34e-02 1.926069e-06 -2.00 3.267549e-05 
-1.72 2.500277e-06 -1.96


                               Pressure convergence table
==
level    dt     cells  dofs      hmax L2               
   H1 Linfty
    6 1.00e-01 10979 20609 1.34e-02 2.855089e-03     - 1.235457e-02   
  -    7.053574e-03     -
    6 2.50e-02 10979 20609 1.34e-02 1.983995e-04 -1.92 2.829075e-03 
-1.06 4.839667e-04 -1.93
    6 6.25e-03 10979 20609 1.34e-02 1.308016e-05 -1.96 2.705119e-03 
-0.03 3.305446e-05 -1.94
    6 1.56e-03 10979 20609 1.34e-02 3.290324e-06 -1.00 2.704622e-03 
-0.00 2.323934e-05 -0.25


The Fig. 4 of the paper presents the log-log plot of the infinity norm 
of the pressure vs time step size using 8 time step sizes. On the 
square domain it shows convergence order of 1.6 on the square domain. 
Whereas on the 

Re: [deal.II] Strong scaling issue on a deal.II based framework running on Skylake nodes.

[Martin Kronbichler]

Dear David,

Without knowing the exact components of deal.II you are using, the first 
places where I would start looking into is whether you use 
multi-threaded blas or multithreading within deal.II. So you could try to do


export DEAL_II_NUM_THREADS=1
export OMP_NUM_THREADS=1

or disable multithreading from the compilation of deal.II (and use 
serial BLAS/LAPACK libraries) and check again. The behavior you're 
describing looks to be a combination of something that sees a good 
speedup in some parts of the solver, but very little to none in other parts.


The second suspicion would be memory bandwidth limitations within the 
node, but even if you are fully memory bound you should see a factor of 
~10-12 of speedup when going from 1 to 48 cores on a node (or a bit less 
if the processor has full turbo frequency turned on and thus clocks 
higher with 1 core loaded than with all 24 cores loaded per socket), 
while you observe much less than that.


Best,
Martin

On 22.01.21 05:52, David Montiel Taboada wrote:

Hello,

I am using the PRISMS-PF framework (which is based on deal.II) on the 
Skylake (skx) nodes (with 48 processors each) of the Stampede2 cluster.


I recently ran a series of strong scaling tests and noticed that the 
intra-node performance (i.e. 1 node, 1-48 processors) scales poorly, 
specifically the solver part. However, once I get past one node, the 
scaling is closer to ideal (taking 1 node as a reference).


Here is the behavior I got (solver part only; in every case I used as 
many MPI threads as processors):


Processors, Nodes, Solver time (s)
1, 1, 821
2 , 1, 608
4 , 1,  525
8, 1, 482
24, 1, 435
48, 1, 427
96, 2, 211
192, 4, 109

Does anyone know what may be the problem?

The code uses  the matrix-free method and requires only the p4est and 
mpi libraries, which I included as dependencies when I did cmake to 
install deal.II.


Here is the line I used
 cmake -DDEAL_II_WITH_MPI=ON -DDEAL_II_WITH_P4EST=ON
 -DCMAKE_INSTALL_PREFIX=$WORK/dealii_install $WORK/dealii-9.2.0

Am I perhaps missing a flag?

By the way, the home nodes (which I used to install deal.II and 
compile my code) are also Skylake, so I would expect my code to have a 
good performance.


I do not observe the same issue elsewhere (e. g., on my local machine 
or on the KNL nodes on Cori).


Any help that might help me figure out this issue is appreciated.

Best,

David




--
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/8bd2837d-c284-4f1e-a194-ad4a56835cb6n%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/b22ad7fa-ccb8-17da-4c5c-daebfe30dd1d%40gmail.com.

Re: [deal.II] Understanding loops in matrix-free

[Martin Kronbichler]

Dear David,

Indeed the two approaches should give the same answer; apart from when 
exactly the cells and faces are scheduled, i.e., possibly a change in 
the order of roundoff errors. Thus, I suspect there is something strange 
going on and could be a bug either in deal.II or in how you set up some 
particular ingredients. Can you share the code or explain a bit what 
exactly is used? We have continuous elements, but which degree? Is it a 
system? How many cells? Are you using threads, i.e., is any of the task 
parallel settings turned on? From what you say I believe you are running 
in serial, right? So it seems a pretty small problem, which we should be 
able to investigate rapidly.


One thing you could do is to look into the "face_range" that you obtain 
when the algorithm calls back into local_apply_boundary_face and compare 
that with the range that you manually construct in your first version? I 
wonder if there are some parts of the loop we are missing or running twice.


Best,
Martin

On 06.02.21 19:18, 'David' via deal.II User Group wrote:


Sorry for messing up the topic. I should be Understanding loops in 
matrix-free. I wanted to insert a figure of the source code rather 
than the google groups formatting and it didn't work for some reason.

David schrieb am Samstag, 6. Februar 2021 um 19:15:07 UTC+1:

Hi there,

I'm currently trying to pack my cell operations into one of the
matrix-free loop functions.
In my first version, I implemented the loops manually by using
(sorry for the odd formatting):
```
    local_apply_cell(*matrix_free,
 system_rhs,
 source_vector,
 std::make_pair(0,
matrix_free->n_cell_batches()));
    local_apply_boundary_face(
  *matrix_free,
  system_rhs,
  source_vector,
std::make_pair(mf_data_reference->n_inner_face_batches(),
mf_data_reference->n_inner_face_batches() +
mf_data_reference->n_boundary_face_batches()));
```
which works as desired. I don't have any operations on internal
faces and ghost exchange/compression doesn't matter for the rest
of the question.

In a second step, I replaced the manually implemented loops by the
matrix-free loop in the following way:
```
    matrix_free->template loop(
  [&](const auto ,
  auto &  dst,
  const auto ,
  const auto _range) {
    local_apply_cell(data, dst, src, cell_range);
  },
  [](const auto &, auto &, const auto &, const auto &) {},
  [&](const auto ,
  auto &  dst,
  const auto ,
  const auto _range) {
    local_apply_boundary_face(data, dst, src, face_range);
  },
  system_rhs,
  source_vector,
  true,
  MatrixFree::DataAccessOnFaces::
    none,
  MatrixFree::DataAccessOnFaces::
    none);
```

However, the system starts diverging after a time (as opposed to
the 'manual' loop), i.e., the result is different.

I should probably add as a comment that I use sorted boundary
condition according to the boundary cell batches so that I apply
the first set of data to the first boundary face batch and the
second one to the second... Though, I checked the ordering
manually and the loop-function seems to work in the same order on
the batches as the manual loop.

I also tried to set the vector zeroing on the destination vector
to 'false' (I use continuous elements) but it doesn't make any
difference. According to my understanding of the loop function,
both implementations should be equivalent. I probably miss
something here.

Maybe anyone else has an idea about the differences.
Thanks in advance and best regards,
David

--
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/620dc4c7-4adf-47f7-8f43-cc9a3dd3e22en%40googlegroups.com 
.


--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 

Re: [deal.II] Incremental build of MG objects

[Martin Kronbichler]

Hi Luca,

Sorry for the delayed answer.


   +--+
   |  |
   |  |
   |Estimate/mark/refine  |
   | ^|
   | +|
   |   Current-level   Current-level  |
   |+^|
   |  Pre smooth|| Post smooth|
   |v+|
   |   Coarser level   Coarser level  |
   |+^|
   |  Pre smooth|| Post smooth|
   |v+|
   |   Coarsest level +--> Coarse Solve   |
   |  |
   +--+

This is an interesting approach.

The issue I’m facing is the following: on each level, I’m setting up the whole 
hierarchy of the MG framework, just to apply one or two V-cycle steps, and the 
whole setup cost is basically comparable to the application of the single 
V-cycle algorithm.

The thing is, once I finished the V-cycle, estimated the error, and transferred 
the solution, the whole hierarchy of MG is destroyed and rebuilt from scratch. 
While this is fine if you are using MG as a preconditioner, and need to call it 
many times, in my case this is the most time consuming part (!).


I can easily imagine that. The question is: Have you measured what 
exactly consumes the time? If you through the whole process again, I 
think the heaviest parts are typically the setup of the MatrixFree 
objects, followed by MGTransferMatrixFree (I assume those are the data 
structures you are using). But the setup of the Triangulation and 
DoFHandler::distribute(_mg)_dofs() should also be quite expensive.




Is there a way to reuse as much as possible of the existing MG objects, i.e., 
detect what levels need to be rebuilt and what can remain the same?


There is a way, but it is not straight-forward and takes a few days of 
work (if you already know where to look). As you suspect, the coarse 
hierarchy depends on the partitioning on the finest level, so the 
objects really need to be refreshed because the parallel distribution 
might change. I think the easiest way would be to hook things up on the 
global coarsening infrastructure that Peter (on CC) has been building 
recently, because there we have already separate triangulations and 
other things. At the same time, the transfer would probably be a bit 
orthogonal to what is done right now, because you would only transfer 
into new unknowns on a few isolated mesh cells.


I think we could definitely identify some strategy to make this work, as 
I see you point of reusing a coarse AMG, for example. How is your time 
plan and resources on this topic?


Best,
Martin

--
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/36a08e9e-2acc-fd1e-5c23-9077f40b457c%40gmail.com.

Re: [deal.II] Tags and Communicators

[Martin Kronbichler]

On 2/19/21 2:38 PM, Timo Heister wrote:

I see additional problems with separate communicators per object:
1. For me it is unclear how operations that involve more than one
object should communicate. For example, a mat-vec has 2 vectors (src,
dst) and a matrix and as such 3 communicators. Of course you can make
an arbitrary choice here, but this doesn't seem clean to me.


That's true, but the point is to avoid different objects communicating 
on the same communicator. Which of the three you mention here is, in 
the end, not important -- in fact, we don't care: Let PETSc or 
Trilinos do what they need to do.


Well, there are parts where we do care, because we have our own parallel 
vectors and parallel matrix-vector products (mostly matrix-free 
functionality). Those cases are linked to the triangulation's 
communicator, which might be the place where duplicating a communicator 
makes sense. One could do it also with a finer granularity inside 
certain collections of vectors, but the crucial ingredient is to make 
data structures aware of the places where objects with different 
communicators come together, like different communicators in the domain 
and range partitioner of a matrix-vector product. It's actually on my 
todo list to use different communicators on coarser levels of the 
multigrid hierarchy, either in terms of different triangulations for 
(global coarsening) or just the levels in the regular h-multigrid with 
local smoothing we have had for a long time.


The important message is that we can't simply go through the library and 
switch over to duplicating communicators for every class where we hand 
in and store a communicator, but, at least for the linear algebra part 
like parallel vectors and MatrixFree, decide on a case-by-case basis so 
that we don't accidentally destroy things.






2. MPI_comm_dup can involve communication and as such can be slow for
large processor counts. Not sure you want to always pay the price for
that. I remember at least one MPI implementation that does an
allreduce (or something worse than that).


That is true too, though one would assume that whatever one ends up 
doing on these communicators is going to be more expensive than the 
duplication.


It depends on what you are doing and on what scale you are running. For 
some MPI implementations, the duplication can be pretty expensive (on 
the order of seconds) once you reach 100k or more ranks, depending on 
how they implement the operation. Imagine what happens if you do an 
MPI_Comm_dup every time you enter a Schur complement approximate inverse 
and request a temporary vector for the local solvers - an operation that 
otherwise scales well down to 1e-3 to 1e-4 seconds. Again, you need a 
case-by-case decision, and everything the below the level of the 
triangulation needs to be audited at least.



It may not be an important issue. Would a compromise be to duplicate 
communicators in the constructors of all main classes in the tutorial, 
as a recommendation of best practices. (This is what ASPECT does, for 
reference.) This doesn't address the issue of what can go wrong when 
using multiple threads from within the main class, but it would at 
least address the issue if someone wanted to run multiple instances of 
the main class in parallel, and may bring the issue to people's mind.


Apart from the disagreement I expressed above, here I completely agree 
with you.


Best,
Martin

--
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/9dce8c01-158f-30a3-f4aa-0a4f70dc5d3e%40gmail.com.

Re: [deal.II] point_value, Real parts, step-29

[Martin Kronbichler]

Dear Hermes,

Did you try VectorTools::point_values, 
https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#a1ad6eceb6cbeaa505baf7f938289bbde 
?


As an alternative, you can also check if the point is in the locally 
owned cell of the respective cell, and run the evaluation only on that 
MPI rank. You can then use an MPI collective, like MPI_Bcast or 
deal.II's version Utilities::MPI::broadcast, to get the solution to all 
MPI ranks (typically with another collective to make all ranks agree on 
who sends the information).


Best,
Martin

This is a variant that allows the evaluation even on remote processes.

On 24.08.21 14:05, Hermes Sampedro wrote:

Dear all,

I am trying to extend my implementation with parallel computing. I 
would like to get the solution in a single point when having a 
distributed system. For example, in step-40 I have the 
/locally_relevant_solution/ at the output_result(...) function.

In a non distributed implementation I used:

Vector<*double*> vecSol (fe.n_components());

VectorTools::point_value(dof_handler, solution,Point<2>(0.2, 0.2),vecSol);


What is the easiest way to do that? Is it possible to check if the 
point is inside the /locally_relevant_solution? /Or would be more 
appropriate to create a copy of the solution in a global vector and 
get the point from there?



Thank you

H


El mié, 11 ago 2021 a las 10:12, Hermes Sampedro 
(mailto:hermesampe...@gmail.com>>) escribió:


Thank you very much, I could make it work.

Regards,
Hermes

El mar, 10 ago 2021 a las 19:46, Jean-Paul Pelteret
(mailto:jppelte...@gmail.com>>) escribió:

Another thing: You probably also need to initialise with the
right number of components. So something like
Vector vecSol (fe.n_components());



On 10. Aug 2021, at 19:40, Jean-Paul Pelteret
mailto:jppelte...@gmail.com>> wrote:

Hi Hermes,

You don’t say what errors you’re seeing, but my guess is that
it now doesn’t compile. This variant of the function (the one
that Daniel linked to) returns void, so you should call it
before outputting the result:

Vector vecSol;
VectorTools::point_value(dof_handler,
solution,Point<2>(0.2, 0.2),vecSol)
std::cout << "Solution at (0.2,0.2): "<< vecSol << std::endl;

This will hopefully get you what you want.

Best,
Jean-Paul



On 10. Aug 2021, at 16:09, Hermes Sampedro
mailto:hermesampe...@gmail.com>>
wrote:

Thank you for your answer. I am not sure if I fully
understand your suggestion. Do you mean something like that:

Vector<*double*> vecSol;
std::cout <<"Solution at (0.2,0.2): "<<
VectorTools::point_value(dof_handler,
solution,Point<2>(0.2,0.2),vecSol)<< std::endl;


I still get some errors. Is there not any way to get for
example the real part of/solution/ easily and use it
directly on the point_value as in step-3?

Thank you
H

El mar, 10 ago 2021 a las 15:49, Daniel Arndt
(mailto:d.arndt.m...@gmail.com>>)
escribió:

Hermes,

Use another overload. The one returning the solution as
a parameter should

work:https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#acd358e9b110ccbf4a7f76796d206b9c7



Best,
Daniel

Am Di., 10. Aug. 2021 um 09:41 Uhr schrieb Hermes
Sampedro mailto:hermesampe...@gmail.com>>:

Dear all,

It is explained in Step-3 how to evaluate the
solution in a point. I am trying to do the same for
Step-29, to evaluate the real and imaginary parts
separately in a single point:

/std::cout << "Solution at (0.2,0.2): "<<
VectorTools::point_value(dof_handler,
solution,Point<2>(0.2, 0.2))<< std::endl;/

I do not have any problems compiling however, an
error occurs when running:

/The violated condition
was:  dof.get_fe(0).n_components() == 1/

What is the proper way to call the real and
imaginary parts of the solution at a particular
point here?


Thank you very much!

H.


--
The deal.II project is located
athttp://www.dealii.org/ 
For mailing list/forum options,
seehttps://groups.google.com/d/forum/dealii?hl=en

---
You received this message because you are subscribed
to the Google 

Re: [deal.II] Post Processor for a mutli-component system

[Martin Kronbichler]

I did not check all code in detail, but what I observe is that you write 
into the joint_solution vector without ghosts. Now this would likely 
result in a run time error (for countinuous elements and in debug mode), 
rather than wrong output. But either way, I would recommend you move the 
"locally_relevant_joint_solution" and the way you initialize it to the 
place where you declare "joint_solution" now, fill the data, and instead 
of calling "compress" you simply call "update_ghost_values".


If that does not help, could you please specify in more detail where the 
output is wrong: Is it just cells near processor interfaces? What kinds 
of elements are you using and what kind of information is the 
postprocessor computing?


Best,
Martin


On 23.09.21 11:40, SALMAN wrote:


Dear all

I am using dealii to solve a multi component system has its own fe, 
dofHandler, and solution vector. I am using a matrix free 
implementation for the entire program.


I am using the Postprocessor class to derive quantities that depend on 
all of the components.


The Postprocessor works fine in serial implementation. However, when 
trying to run using MPI the output is producing wrong results.


I created a procedure to generate a joint solution vector and a joint 
handler following tutorial 32. However, this tutorial uses Trillinos 
vectors and I am using dealii’s own distributed vector and I am 
suspecting that this might be problem.


Hence, my question is how to properly generate a joint solution vector 
while using dealii's distributed::vector class.



I have attached a text file that contains the procedure that was used 
to generate a the joint vector and I would appreciate if anyone would 
point out the mistakes that I made



Thank you

​
--
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/c4beae48-c449-4b03-85df-1406bebc11dbn%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/71317c8d-a236-ccf5-3b03-0ebb3a1d95da%40gmail.com.

Re: [deal.II] Re: Interlacing update ghost and assemble

[Martin Kronbichler]

Dear Praveen,

In addition to what Bruno said, I would just add that I would benchmark 
the respective costs carefully before putting it into the code. As you 
observed, this requires some additional steps to figure out what work 
can be reasonably overlapped. We have this in the matrix-free framework 
where enough interest has been available, but not for the rest (where we 
often did not have overly performance-sensitive applications).


But even if you are extremely performance-sensitive, I would not 
immediately subscribe to a big benefit. My observation has been that it 
may bring a few percent, but not big factors, and much less what 
"general wisdom" seems to indicate (and you can hardly find papers that 
quantify the benefit in a reproducible way). I guess networks have 
improved quite a bit as to when this wisdom was shaped. Done naively and 
especially on CPU systems it can even cause slowdowns: When you have the 
bulk of MPI communication within a node (rather than across nodes), it 
is typically the memory bandwidth that is limiting. You would hence just 
shift that to another stage in your code, at the expense of fooling the 
prefetchers by non-standard loop through cells. I found that this could 
cost more than it helps in some contexts. Things are different when you 
have lots of inter-node communication (say over Infiniband) and you can 
do more useful work, or when you have GPUs where the computations are 
quicker in general. But the biggest success stories I have seen are 
really complicated to implement, needing separate cores for the MPI 
communication as opposed to other "worker threads" and other tricks. 
Again, I am not disputing it can be worth it, but I would first look 
whether that is one of the low-hanging fruits or just disappears in the 
noise.


Best,
Martin

On 01.10.21 14:46, Bruno Turcksin wrote:

Praveen,

We do something like that in CUDA MatrixFree. It is slightly more 
complicated because we need to update the ghost values on both the 
source vector and the destination vector. The idea is to first loop 
over the mesh and store all the vertices that are ghosted then you 
loop over the active cells, you just need to check that none of the 
vertices is ghosted. The relevant code is here: 
https://github.com/dealii/dealii/blob/master/include/deal.II/matrix_free/cuda_matrix_free.templates.h#L1012-L1057


Best,

Bruno

On Thursday, September 30, 2021 at 11:48:13 PM UTC-4 Praveen C wrote:

Dear all

I use la::d::Vector in a continuous nodal FE code. I want to mix
update ghost and assemble so that I dont wait for update ghost to
finish

solution.update_ghost_values_start()
Assemble on all locally owned cells that do not need ghost values
solution.update_ghost_values_finish()
Assemble on all locally owned cells that need ghost values

I use WorkStream for assembly.

I have to create a FilteredIterator that identifies all cells that
do not have ghost dof values.

One way is to check all dofs on a cell to see if all of them are
locally owned.

Is there any other way to create such a filter ?

I have to use two WorkStreams here, the second one will have less
work to do. Is there a way to avoid using two WorkStreams ?

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/8e65c3ba-be99-497d-ad62-1ea226cace90n%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/7a00a20e-a680-517c-6827-ac2a5c019199%40gmail.com.

Re: [deal.II] Are there known issue between Trilinos 13.0.1 and deal.II?

[Martin Kronbichler]

Dear Bruno,

This is a regular error and not related to the actual problem I believe. 
This is simply the message that we use to detect what SIMD extensions we 
can enable, and in your case the compilation does not support AVX. But I 
think further down in the CMakeFiles/CMakeError.log file you should find 
the culprit.


As a side remark, the AMD Rome architecture does indeed support AVX2, as 
long as you pass "-mavx2 -mfma" or "-march=znver2" or simply 
"-march=native" to the compile flags ("CMAKE_CXX_FLAGS"). However, AMD 
does not (yet) support AVX-512, which might be what you're referring to. 
In other words, the compiler needs to be told what kind of architecture 
extensions are supported on the target CPU. I most often use 
"-march=native" unless I'm cross-compiling on a different CPU compared 
to where the code is eventually executed.


Best,
Martin

On 29.10.21 12:47, blais...@gmail.com wrote:

Wolfgang, you are a genius.
The exact error I am getting is related to AVX instructions not being 
available. I think this is normal because this is an AMD Rome based 
cluster and not an intel one.


Run Build 
Command(s):/cvmfs/soft.computecanada.ca/gentoo/2020/usr/bin/gmake 
cmTC_be568/fast && 
/cvmfs/soft.computecanada.ca/gentoo/2020/usr/bin/gmake -f 
CMakeFiles/cmTC_be568.dir/build.make CMakeFiles/cmTC_be568.dir/build
gmake[1]: Entering directory 
'/home/blaisbru/dealii/build/CMakeFiles/CMakeTmp'

Building CXX object CMakeFiles/cmTC_be568.dir/src.cxx.o
/cvmfs/soft.computecanada.ca/easybuild/software/2020/Core/gcccore/9.3.0/bin/c++ 
-DDEAL_II_HAVE_AVX   -o CMakeFiles/cmTC_be568.dir/src.cxx.o -c 
/home/blaisbru/dealii/build/CMakeFiles/CMakeTmp/src.cxx
/home/blaisbru/dealii/build/CMakeFiles/CMakeTmp/src.cxx:3:6: error: 
#error "__AVX__ flag not set, no support for AVX"

    3 | #error "__AVX__ flag not set, no support for AVX"
  |  ^
/home/blaisbru/dealii/build/CMakeFiles/CMakeTmp/src.cxx: In function 
‘int main()’:
/home/blaisbru/dealii/build/CMakeFiles/CMakeTmp/src.cxx:35:9: warning: 
AVX vector return without AVX enabled changes the ABI [-Wpsabi]

   35 |   b = _mm256_set1_pd (static_cast(2.25));
  | ~~^
gmake[1]: *** [CMakeFiles/cmTC_be568.dir/build.make:66: 
CMakeFiles/cmTC_be568.dir/src.cxx.o] Error 1
gmake[1]: Leaving directory 
'/home/blaisbru/dealii/build/CMakeFiles/CMakeTmp'

gmake: *** [Makefile:121: cmTC_be568/fast] Error 2

Hopefully I find a solution from there :)!


On Thursday, October 28, 2021 at 11:20:42 p.m. UTC-4 Wolfgang Bangerth 
wrote:


On 10/28/21 8:02 PM, blais...@gmail.com wrote:
>
> However, whenever I run CMAKE using deal.II AND enabling
Trilinos, I get the
> following error:
> Configuration error: Cannot compile a test program with the
final set of
>     compiler and linker flags:
>     CXX flags (DEBUG): -pedantic -fPIC -Wall -Wextra
-Wmissing-braces
> -Woverloaded-virtual -Wpointer-arith -Wsign-compare
-Wsuggest-override
> -Wswitch -Wsynth -Wwrite-strings -Wno-placement-new
> -Wno-deprecated-declarations -Wno-literal-suffix -Wno-psabi
> -Wno-class-memaccess -fopenmp-simd -Wno-parentheses
-Wno-unused-local-typedefs
> -O0 -ggdb -Wa,--compress-debug-sections

Somewhere, in some file that is probably located under CMakeFiles
in your
build directory, it will provide you with whatever error messages
cmake got
when it tried to compile a test program with these flags. If you
know what
that error is, you're probably 75% there towards identifying what
the real
underlying problem is.

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/fb910de8-53f0-4eca-90a7-14a54e3069a1n%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 

Re: [deal.II] shape function gradient at arbitrary points in the element

[Martin Kronbichler]

Dear Simon,

We have the class FEPointEvalation, 
https://dealii.org/developer/doxygen/deal.II/classFEPointEvaluation.html 
, which implements the operation of evaluating a solution at arbitrary 
points (as handed in as an array to points in unit coordinates). To do 
this for positions in real coordinates, you typically need to invert the 
mapping, Mapping::transform_real_to_unit_cell or 
Mapping::transform_points_real_to_unit_cell, see 
https://dealii.org/developer/doxygen/deal.II/classMapping.html , always 
assuming that you have located the correct points. Most tasks can be 
reduced to this setup.


If your intent is not just the solution interpolation but the value of 
shape functions, e.g. because you want to assemble a local matrix, 
FEPointEvaluation is not directly providing that info, so in that case 
the recommendation would be to use FEValues, despite that being terribly 
inefficient. We could possibly integrate such functionality into 
FEPointEvaluation if you have interest in that direction, so let us know 
what your needs are.


Best,
Martin

On 05.07.22 15:56, Simon wrote:

Dear all:

I have to compute the gradient with respect to real space co-ordinates 
of the shape functions belonging to a FE_Q element at arbitrary points 
in the element, that is, not only at quadrature points.


What is the way to do this in dealii?
FEValues provides shape function values, gradients,... only at 
quadrature points.


Thank you,
Simon
--
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/16afc8fa-bcfe-4a81-9ade-040c18f85596n%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/34843cd0-22a6-d4b5-7cb8-f530715a6108%40gmail.com.

Re: [deal.II] Re: Run time analysis for step-75 with matrix-free or matrix-based method

[Martin Kronbichler]

Dear Wayne,

I am a bit surprised by your numbers and find them rather high, at least 
with the chosen problem sizes. I would expect the matrix-free solver to 
run in less than a second for 111,000 unknowns on typical computers, not 
almost 10 seconds. I need to honestly say that I do not have a good 
explanation at this point. I did not write this tutorial program, but I 
know more or less what should happen. Let me ask a basic question first: 
Did you record the timings with release mode? The numbers would make 
more sense if they are based on the debug mode.


Best,
Martin


On 19.10.22 12:08, 'yy.wayne' via deal.II User Group wrote:

Thanks for your reply Peter,

The matrix-free run is basic same as in step-75 except I substitute 
coarse grid solver. For fe_degree=6 without GMG and fe_degree in each 
level decrease by 1 for pMG, the solve_system() function runtime is 
24.1s. It's decomposed to *MatrixFree MG operators 
construction*(1.36s), MatrixFree MG transfers(2.73s),  KLU coarse grid 
solver(5.7s), *setting smoother_data and compute_inverse_diagonal for 
level matrices*(3.4s) CG iteration(9.8s).


The two bold texts cost a lot more(133s and 62s, respectively) in 
matrix-based multigrid case. I noticed just as in step-16, the finest 
level matrix is assembled twice(one for system_matrix and one for 
mg_matrices[maxlevel]) so assembling time cost more.


Best,
Wayne

在2022年10月19日星期三 UTC+8 17:10:27 写道：

Hi Wayne,

your numbers make totally sense. Don't forget that you are running
for high order: degree=6! The number of non-zeroes per
element-stiffness matrix is ((degree + 1)^dim)^2 and the cost of
computing the element stiffness matrix is even ((degree +
1)^dim)^3 if I am not mistaken (3 nested loop: i, j and q). Higher
orders are definitely made for matrix-free algorithms!

Out of curiosity: how large is the setup cost of MG in the case of
the matrix-free run? As a comment: don't be surprised that the
setup costs are relatively high compared to the solution process:
you are probably setting up a new Triangulation-, DoFHander-,
MatrixFree-, ... -object per level. In many simulations, you can
reuse these objects, since you don't perform AMR every time step.

Peter

On Wednesday, 19 October 2022 at 10:38:34 UTC+2 yy.wayne wrote:

Hello everyone,

I modified step-75 a little bit and try to test it's runtime.
However the result is kind of inexplainable from my point of
view, especially on *disproportionate assemble time and solve
time*. Here are some changes:
1. a matrix-based version of step75 is contructed to compare
with matrix-free one.
2. no mesh refinement and no GMG, and fe_degree is constant
across all cells within every cycle. Fe_degree adds one after
each cycle. I make this setting to compare runtime due to
fe_degree.
3. a direct solver on coareset grid. I think it won't affect
runtime since coarest grid never change

For final cycle it has fe_degree=6 and DoFs=111,361.
For matrix-based method, overall runtime is 301s where setup
system(84s) and solve system(214s) take up most. In step-75
solve system actually did both multigrid matrices assembling,
smoother construction, and CG solving. Runtime of this case is
shown:
matrix-based.png
On each level I print time assembling level matrix. *The solve
system is mostly decomposed to MG matrices
assembling(83.9+33.6+...=133s), smoother set up(65s), coarse
grid solve(6s) and CG solve(2.56).* My doubt is why actual CG
solve only takes 2.56 out of 301 seconds for this problem? The
time spent on assembling and smoother construction account too
much that they seems a burden.

For matrix-free method however, runtime is much smaller
without assembling matrices. Besides, CG solve cost more
because of more computation required by matrix-free I guess.
But *smoother construction time reduces significantly* as well
is out of my expectation.
matrix-free.png

Matrix-free framework saves assembling time but it seems too
efficient to be real. The text in bold are my main confusion.
May someone share some experience on matrix-free and multigrid
methods' time consumption？

Best,
Wayne

--
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/58dae75e-644e-49d8-bee2-e212c5184e1cn%40googlegroups.com