Re: [deal.II] Solving time dependent heat equation with MPI (PETsc)

[Wolfgang Bangerth]

On 10/16/2017 01:53 PM, Mark Ma wrote:


I think this problem lies in the time updating of solution using 
old_solution, since the mass_matrix and laplace_matrix have already 
eliminated the constraint node, /*mass_matrix_T.vmult 
(system_rhs, old_solution_T_cal*//*);*/ is no longer valid for this 
case.


Yes, this sounds reasonable from your description. Have you looked at 
steps 24, 25, 26 to see how it is done there? I could imagine that you 
need to think about which degrees of freedom you want eliminated/not 
eliminated in the matrices you multiply with. To this end, write out the 
weak form of the problem you need to solve in each time step, and think 
specifically about the boundary values of the trial and test functions 
and whether they should be part of the matrix you use on the right hand 
side or not. It is possible that you need to use a different matrix for 
the lhs and rhs operations.


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.

Re: [deal.II] Copy blocks of BlockSparseMatrix

[Alex Jarauta]

Hi Daniel,

thanks for your kind reply. I apologize for not giving more details about 
my code.

In the file where the global matrix (m33) and its sparsity pattern is 
created, I also include the declaration and initialization of the new 
matrix (m22). After that, the new matrix is of size 2x2 blocks. Then, I 
simply copy the sparsity pattern of each block of matrix m33 to the blocks 
of matrix m22 as:

for (int i = 0; i < dim; i++)
{
for (int j = 0; j < dim; j++)
{
m22.block(i,j).reinit(m33.block(i+1,j+1).get_sparsity_pattern()); 
m22.block(i,j).copy_from(m33.block(i+1,j+1));
}
}

The error appears on the next step, when I try to initialize the 
preconditioner with matrix m22:

SparseDirectUMFPACK preconditioner;
preconditioner.initialize (A_mat, SparseDirectUMFPACK::AdditionalData());

The exact error message that I get in debug mode is the following:


An error occurred in line <272> of file 

 
in function
void dealii::SparseDirectUMFPACK::factorize(const Matrix&) [with Matrix 
= dealii::BlockSparseMatrix]
The violated condition was: 
row_pointers[i] == Ap[i+1]
The name and call sequence of the exception was:
ExcInternalError()
Additional Information: 
This exception -- which is used in many places in the library -- usually 
indicates that some condition which the author of the code thought must be 
satisfied at a certain point in an algorithm, is not fulfilled. An example 
would be that the first part of an algorithm sorts elements of an array in 
ascending order, and a second part of the algorithm later encounters an an 
element that is not larger than the previous one.

There is usually not very much you can do if you encounter such an 
exception since it indicates an error in deal.II, not in your own program. 
Try to come up with the smallest possible program that still demonstrates 
the error and contact the deal.II mailing lists with it to obtain help.


I cannot understand why the initialization of the blocks of matrix m22 is 
not working properly..

Thanks for your help.

Cheers,


Alex

On Monday, October 16, 2017 at 9:03:05 AM UTC-6, Daniel Arndt wrote:
>
> Alex,
>
> Hi all,
>>
>> I am having a similar problem. I want to take a 3x3 BlockSparseMatrix m33 
>> and copy four of its blocks into a 2x2 BlockSparseMatrix m22. For instance:
>>
>> m33 = [ A B C
>>  D E F
>>  G H I ]
>>
>> From this matrix, I would like to obtain the following:
>>
>> m22 = [ E F
>>  H I ]
>>
>> Timo, I have tried the solutions that you suggest, but none of them work, 
>> I keep getting a segmentation fault error. Could you please give me an 
>> alternative option?
>>
> What exactly are you doing and what is the exact error message you are 
> getting in debug mode?
> Do the sparsity patterns already match before copying?
>
> 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.

Re: [deal.II] Need for simple applied examples

[Mark Ma]

Hi Konstantin,

Very very long time not contact with you, how are you doing? I am sorry I 
really have not too much time that could spend on Deal.II project, the 
Maxwell equation's problem and so on. 

what you propose here I think is right, at least in the view of very 
beginner of deal.II like me, the documentation in deal.ii is not hard to 
start with but really hard to understand and implement more advanced 
techniques like Adaptive mesh, hp refinement, PETsc / Trilinos. Since these 
techniques are introduced in different problems, very beginners may be 
expertise in optics but have to learn these techniques from tutorials of 
fludic-dynamics. Although from mathematics point of view, this is simple. I 
believe if there are some kind of simple examples implementing these 
techniques gradually would greatly decrease the learning and thinking time 
of beginners.

By the way, right now I am calculating temperature in MPI and unfortunately 
stucked. 

Best,
Mark
--
Laboratoire Hubert Curien, UMR CNRS 5516,
Bâtiment F 18 Rue du Professeur Benoît Lauras
42000 Saint-Etienne
FRANCE

在 2015年9月7日星期一 UTC+2下午11:26:55，Konstantin Ladutenko写道：
>
> Hi,
>
> Sorry, it was hard for me to find enough time to get really started with 
> dell. Few things I was able to learn this summer you can find in this group 
> and in updated step-6.
>
> The next possible step for me is to solve Mie problem using deal.ii and 
> compare against analytic solution. It something similar to H. Fahs 
> self-written code you can find at http://hfahs.free.fr/software.html (and 
> related bunch of publications http://hfahs.free.fr/pubs.html ). For sure 
> I will present it here as soon as I have something to solve the problem of 
> EM scattering correctly (and probably after that it can be added some 
> nonlinearity for plasmonic nanoparticles). 
>
> I also found a book of Peter Monk "FEM for Maxwell equations" that is 
> rather hard for me to read, however, I hope, it will give me answers to all 
> my basic questions. There is also "Introduction to the Finite Element 
> Method in Electromagnetics" by
> Polycarpou which is much easier to read from my "applied" point of view.
>
> Best regards,
> Kostya
>
>
> пн, 7 сент. 2015 г. в 21:21, JT Hu :
>
>> Hi Kostya,
>>
>> It has been a year since you posted this.  Do you have any updates on 
>> developing Deal II methods for plasmonics/photonics calculations to share? 
>> If you have published any related work discuss the methods and coding, I 
>> will also be glad to read and cite!
>>
>> Best,
>> Jingtian
>>
>>
>> On Tuesday, September 16, 2014 at 4:44:42 AM UTC-5, Konstantin Ladutenko 
>> wrote:
>>
>>> Hi Wolfgang,
>>>
>>>  Thank you for your reply! Deal.II documentation is great, I really like 
>>> your lectures, I will recommend them to my students.
>>>
>>> 1)  My goal is to correctly solve my problem as fast as possible, time 
>>> needed to learn new tool included.
>>>
>>> 2) My choice is open source. We have few Comsol licenses usually not 
>>> available for students (taken out for 'great and terrible' top science 
>>> projects), this also leads to some troubles than you apply to use external 
>>> cluster for computations. There is a number of other open source FEM codes 
>>> (Elmer, GetFEM++, Moose, etc.), however as to me deal.ii seems to be a 
>>> mature tool which I am learning with hope to use it many years after that.
>>>
>>> 3) I suppose that almost any project can benefit from large user base. 
>>> It seems that deal.ii can solve many problems in different research areas. 
>>> However average researcher with need of calculations will prefer to use 
>>> matlab (or octave, scilab), may be some python codes (e.g. with help of 
>>> sage) and sometimes some ancient Fortran codes - too large and awful to be 
>>> maintained. GUI programs (Comsol, Ansys etc) are also an option.
>>>
>>> It is a good way to start from basic elements of deal.ii, as it is done 
>>> in stepXX tutorials and to go to some real world problems step by step. 
>>> However you will need few weeks before your can start with your problem. 
>>> Many possible users of deal.ii do not have so much time to be spend without 
>>> any practical results. May be I am wrong, but a number of simple applied 
>>> examples can be used as a kick start entry point to deal.ii. You can use 
>>> them to copy-paste as building bricks for your exact (may be simplified) 
>>> problem. It will not be solved in the best way deal.ii can provide,  but it 
>>> will motivate to learn more details about how it should be done in deal.ii 
>>> way.
>>>
>>> 4) My current research area of interest is metamaterials (with nonlinear 
>>> effects coming soon), so Maxwell equations are one the first place. I am 
>>> mostly interested in transient problems, however after 20 video lectures I 
>>> am still not sure how i should treat them in deal.ii). This is my entry 
>>> list of examples for electromagnetism (all of them can be 

Re: [deal.II] Solving time dependent heat equation with MPI (PETsc)

[Mark Ma]

source code are attached.

在 2017年10月16日星期一 UTC+2下午9:14:13，Wolfgang Bangerth写道：
>
> On 10/16/2017 09:23 AM, Mark Ma wrote: 
> > 
> > It almost looks to me like you're applying both Dirichlet values (via 
> > the ConstraintMatrix) and Neumann boundary values (via a boundary 
> > integral). But then I haven't taken a look at the code, so I can't 
> > really say for sure. 
> > 
> > 
> > I think I only applied Dirichlet values via ConstraintMatrix. There 
> > is no boundary integrals, so I believe Neumann BC is not 
> > implemented. 
>
> My guess came from the fact that a Neumann boundary condition would act 
> like a boundary heat source. It looks like you have a heat source in the 
> last layer of cells around the boundary. 
>
> I don't know whether that is true, but it's a place I'd try to 
> investigate in debugging. 
>
>
> > But may be you are right since for most of the case, heat equations 
> > are solved assuming an adiabatic process, that is dT/dx = 0 at 
> > boundaries. This zero values is automatically satisfied. 
>
> That's a different question. Whether one may *physically* want a 
> different equation is besides the point -- you have chosen one equation 
> you'd like to solve, but the numerical solution is wrong. It's 
> worthwhile trying to understand why that is before thinking about 
> *other* equations. 
>
>
> > I think the problem comes from the way of setting Boundary values 
> > using constraint method, which in dealii it makes some kind of 
> > constraints at vertex of boundaries, i.e. x22 = x1/2 + x2/2, 
>
> Hm, that would be the constraint of a hanging node. Do you have hanging 
> nodes in your problem? If you do, may I suggest to simplify the problem 
> and see if you have the same problem on uniform meshes? 
>
> 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.
/* -
 *
 * Copyright (C) 1999 - 2016 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.
 *
 * -

 *
 * Author: Wolfgang Bangerth, University of Heidelberg, 1999
 */


// @sect3{Include files}

// The first few (many?) include files have already been used in the previous
// example, so we will not explain their meaning here again.
//#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 
#include 

// This is new, however: in the previous example we got some unwanted output
// from the linear solvers. If we want to suppress it, we have to include this
// file and add a single line somewhere to the program (see the main()
// function below for that):
#include 

#include 

#include 

// The final step, as in previous programs, is to import all the deal.II class
// and function names into the global namespace:
using namespace dealii;


// @sect3{The Step4 class template}

// This is again the same Step4 class as in the previous
// example. The only difference is that we have now declared it as a class
// with a template parameter, and the template parameter is of course the
// spatial dimension in which we would like to solve the Laplace equation. Of
// course, several of the member variables depend on this dimension as well,
// in particular the Triangulation class, which has to represent
// quadrilaterals or hexahedra, respectively. Apart from this, everything is
// as before.
template 
class Step4
{
public:
  Step4 ();
  void run ();

private:
  void make_grid ();
  void setup_system();
  void assemble_system ();

  void assemble_rhs_T(double rhs_time);
  void solve_T ();
  void solve_T_run ();
  void output_results (int output_num) const;
  
  MPI_Comm 

Re: [deal.II] Solving time dependent heat equation with MPI (PETsc)

[Mark Ma]

> On 10/16/2017 09:23 AM, Mark Ma wrote: 
> > 
> > It almost looks to me like you're applying both Dirichlet values (via 
> > the ConstraintMatrix) and Neumann boundary values (via a boundary 
> > integral). But then I haven't taken a look at the code, so I can't 
> > really say for sure. 
> > 
> > 
> > I think I only applied Dirichlet values via ConstraintMatrix. There 
> > is no boundary integrals, so I believe Neumann BC is not 
> > implemented. 
>
> My guess came from the fact that a Neumann boundary condition would act 
> like a boundary heat source. It looks like you have a heat source in the 
> last layer of cells around the boundary. 
>
> I don't know whether that is true, but it's a place I'd try to 
> investigate in debugging. 
>
>
Yes, it looks like so. But then I comment the time updating code of rhs 
source by the following code, this boundary heat still appears(see figs 
below)
 //
  //time dependent

//assign right hand side
mass_matrix_T.vmult (system_rhs, old_solution_T_cal);
laplace_matrix_T.vmult (tmp, old_solution_T_cal);
system_rhs.add (-time_step * (1-theta), tmp);
/*
assemble_rhs_T (time);
forcing_terms = dynamic_rhs_T;
forcing_terms *= time_step * theta;

assemble_rhs_T (time - time_step); 
tmp = dynamic_rhs_T;
forcing_terms.add (time_step*(1-theta),tmp);
system_rhs.add (1,forcing_terms);

*/

//assign system matrix
system_matrix.copy_from (mass_matrix_T);
system_matrix.add (time_step * theta, laplace_matrix_T);

I think this problem lies in the time updating of solution using 
old_solution, since the mass_matrix and laplace_matrix have already 
eliminated the constraint node, 
*mass_matrix_T.vmult (system_rhs, old_solution_T_cal**);* is no longer 
valid for this case. This may result into ZEROs at all constrained nodes, 
so at boundaries these zero values then propagate into center, as can be 
seen from the figure below. For time updating of solution using 
old_solution, maybe tutorial-32 could give some ideas. However to 
completely understand step-32 may spend some more time for me. 






 

>
> > But may be you are right since for most of the case, heat equations 
> > are solved assuming an adiabatic process, that is dT/dx = 0 at 
> > boundaries. This zero values is automatically satisfied. 
>
> That's a different question. Whether one may *physically* want a 
> different equation is besides the point -- you have chosen one equation 
> you'd like to solve, but the numerical solution is wrong. It's 
> worthwhile trying to understand why that is before thinking about 
> *other* equations. 
>
>
> > I think the problem comes from the way of setting Boundary values 
> > using constraint method, which in dealii it makes some kind of 
> > constraints at vertex of boundaries, i.e. x22 = x1/2 + x2/2, 
>
> Hm, that would be the constraint of a hanging node. Do you have hanging 
> nodes in your problem? If you do, may I suggest to simplify the problem 
> and see if you have the same problem on uniform meshes? 
>
>
there is no handing node in my code since I used a uniformed mesh. And also 
I tried this, results show no difference.

I am really appreciate your patience to look for this question during your 
busy schedule.

Thanks and have a nice day,
Mark
-
Laboratoire Hubert Curien, UMR CNRS 5516,
Bâtiment F 18 Rue du Professeur Benoît Lauras
42000 Saint-Etienne
FRANCE

-- 
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: Error when applying initial values to MPI::Vector in multiple dimensions

[Denis Davydov]

interpoalate just evaluates the function at support points of FE basis 
(assuming that you have one with support points)
and sets those values to DoFs,
whereas project (as the name implies) does L2 projection. Thus as others have 
mentioned you are solving
Mx = U 
where M is the mass matrix.

Regards,
Denis.

> On 16 Oct 2017, at 21:44, 'Maxi Miller' via deal.II User Group 
>  wrote:
> 
> Yes, .interpolate() works fine. But what is the difference between 
> interpolate() and project()?
> 
> Am Montag, 16. Oktober 2017 08:54:02 UTC+2 schrieb Denis Davydov:
> Or you may want to use interpolate if this is enough
> https://www.dealii.org/developer/doxygen/deal.II/namespaceVectorTools.html#a05db6c8cebf924b417dd92f525efe3db
>  
> 
> 
> Regards,
> Denis
> 
> On Friday, October 13, 2017 at 11:05:01 AM UTC+2, Maxi Miller wrote:
> I try to apply initial values to a vector defined as 
> LinearAlgebraTrilinos::MPI::Vector using 
> VectorTools::project (dof_handler, hanging_node_constraints,
>  QGauss(fe.degree+1),
>  InitialValues(),
>  local_solution);
> 
> 
> When initializing the variable fe (as FESystem) with one or two 
> components, it works fine. For more than two components I get the error
>  
> An error occurred in line <1366> of file 
> <~/Downloads/dealii/include/deal.II/numerics/vector_tools.templates.h> in 
> function 
> void dealii::VectorTools::{anonymous}::project(const 
> dealii::Mapping&, const dealii::DoFHandler&, const 
> dealii::ConstraintMatrix&, const dealii::Quadrature&, const 
> dealii::Function&, VectorType&, bool, 
> const dealii
> ::Quadrature<(dim - 1)>&, bool) [with VectorType = 
> dealii::TrilinosWrappers::MPI::Vector; int dim = 2; typename 
> VectorType::value_type = double] 
> The violated condition was:  
> (dynamic_cast* > 
> (&(dof.get_triangulation()))==nullptr) 
> Additional information:  
> You are trying to use functionality in deal.II that is currently not 
> implemented. In many cases, this indicates that there simply didn't appear 
> much of a need for it, or that the author of the original code did not have 
> the time to implement a particular case. If you
>  hit this exception, it is therefore worth the time to look into the code to 
> find out whether you may be able to implement the missing functionality. If 
> you do, please consider providing a patch to the deal.II development sources 
> (see the deal.II website on how to contri
> bute). 
>  
> Stacktrace: 
> --- 
> #0  /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre:  
> #1  /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre: void 
> dealii::VectorTools::project<2, dealii::TrilinosWrappers::MPI::Vector, 
> 2>(dealii::DoFHandler<2, 2> const&, dealii::ConstraintMatrix const&, 
> dealii::Quadrature<2> const&, dealii::Function<2, 
> dealii::TrilinosWrappers::MPI
> ::Vector::value_type> const&, dealii::TrilinosWrappers::MPI::Vector&, bool, 
> dealii::Quadrature<(2)-(1)> const&, bool) 
> #2  ./main: Step15::MinimalSurfaceProblem<2>::run() 
> #3  ./main: main 
>  
>  
> [linux-lb8c:15830] *** Process received signal *** 
> [linux-lb8c:15830] Signal: Aborted (6) 
> [linux-lb8c:15830] Signal code:  (-6) 
> [linux-lb8c:15830] [ 0] /lib64/libpthread.so.0(+0x12270)[0x7f294a477270] 
> [linux-lb8c:15830] [ 1] /lib64/libc.so.6(gsignal+0x110)[0x7f2946c1f0d0] 
> [linux-lb8c:15830] [ 2] /lib64/libc.so.6(abort+0x151)[0x7f2946c206b1] 
> [linux-lb8c:15830] [ 3] 
> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(+0x6b9e5d1)[0x7f295b49e5d1] 
> [linux-lb8c:15830] [ 4] 
> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals5abortERKNS_13ExceptionBaseE+0x1a)[0x7f295b49edaf]
>  
> [linux-lb8c:15830] [ 5] 
> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals11issue_errorINS_18StandardExceptions17ExcNotImplementedEEEvNS1_17ExceptionHandlingEPKciS7_S7_S7_T_+0x98)[0x7f2957373ea1]
>  
> [linux-lb8c:15830] [ 6] 
> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(+0x3f38e23)[0x7f2958838e23] 
> [linux-lb8c:15830] [ 7] 
> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii11VectorTools7projectILi2ENS_16TrilinosWrappers3MPI6VectorELi2EEEvRKNS_10DoFHandlerIXT_EXT1_EEERKNS_16ConstraintMatrixERKNS_10QuadratureIXT_EEERKNS_8FunctionIXT1_ENT0_10value_typeEEERSH_bRKNSC_IX
> miT_Li1b+0x2f)[0x7f295894906e] 
> [linux-lb8c:15830] [ 8] 
> ./main(_ZN6Step1521MinimalSurfaceProblemILi2EE3runEv+0xc08)[0x420d08] 
> [linux-lb8c:15830] [ 9] ./main(main+0x3c)[0x414ad0] 
> [linux-lb8c:15830] [10] 
> /lib64/libc.so.6(__libc_start_main+0xea)[0x7f2946c09f4a] 
> [linux-lb8c:15830] [11] ./main(_start+0x2a)[0x41477a] 
> [linux-lb8c:15830] *** End of error message *** 
> Abgebrochen (Speicherabzug geschrieben)
> 
> when running in debug mode. 

[deal.II] Re: Error when applying initial values to MPI::Vector in multiple dimensions

['Maxi Miller' via deal.II User Group]

Yes, .interpolate() works fine. But what is the difference between 
interpolate() and project()?

Am Montag, 16. Oktober 2017 08:54:02 UTC+2 schrieb Denis Davydov:
>
> Or you may want to use interpolate if this is enough
>
> https://www.dealii.org/developer/doxygen/deal.II/namespaceVectorTools.html#a05db6c8cebf924b417dd92f525efe3db
>
> Regards,
> Denis
>
> On Friday, October 13, 2017 at 11:05:01 AM UTC+2, Maxi Miller wrote:
>>
>> I try to apply initial values to a vector defined as 
>> *LinearAlgebraTrilinos::MPI::Vector 
>> *using 
>> VectorTools::project (dof_handler, hanging_node_constraints,
>>  QGauss(fe.degree+1),
>>  InitialValues(),
>>  local_solution);
>>
>>
>>
>> When initializing the variable fe (as FESystem) with one or two 
>> components, it works fine. For more than two components I get the error
>>
>>  
>> An error occurred in line <1366> of file 
>> <~/Downloads/dealii/include/deal.II/numerics/vector_tools.templates.h> in 
>> function 
>> void dealii::VectorTools::{anonymous}::project(const 
>> dealii::Mapping&, const dealii::DoFHandler&, const 
>> dealii::ConstraintMatrix&, const dealii::Quadrature&, const 
>> dealii::Function&, VectorType&, bool, 
>> const dealii
>> ::Quadrature<(dim - 1)>&, bool) [with VectorType = 
>> dealii::TrilinosWrappers::MPI::Vector; int dim = 2; typename 
>> VectorType::value_type = double] 
>> The violated condition was:  
>> (dynamic_cast* > 
>> (&(dof.get_triangulation()))==nullptr) 
>> Additional information:  
>> You are trying to use functionality in deal.II that is currently not 
>> implemented. In many cases, this indicates that there simply didn't appear 
>> much of a need for it, or that the author of the original code did not have 
>> the time to implement a particular case. If you
>>  hit this exception, it is therefore worth the time to look into the code to 
>> find out whether you may be able to implement the missing functionality. If 
>> you do, please consider providing a patch to the deal.II development sources 
>> (see the deal.II website on how to contri
>> bute). 
>>  
>> Stacktrace: 
>> --- 
>> #0  /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre:  
>> #1  /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre: void 
>> dealii::VectorTools::project<2, dealii::TrilinosWrappers::MPI::Vector, 
>> 2>(dealii::DoFHandler<2, 2> const&, dealii::ConstraintMatrix const&, 
>> dealii::Quadrature<2> const&, dealii::Function<2, 
>> dealii::TrilinosWrappers::MPI
>> ::Vector::value_type> const&, dealii::TrilinosWrappers::MPI::Vector&, bool, 
>> dealii::Quadrature<(2)-(1)> const&, bool) 
>> #2  ./main: Step15::MinimalSurfaceProblem<2>::run() 
>> #3  ./main: main 
>>  
>>  
>> [linux-lb8c:15830] *** Process received signal *** 
>> [linux-lb8c:15830] Signal: Aborted (6) 
>> [linux-lb8c:15830] Signal code:  (-6) 
>> [linux-lb8c:15830] [ 0] /lib64/libpthread.so.0(+0x12270)[0x7f294a477270] 
>> [linux-lb8c:15830] [ 1] /lib64/libc.so.6(gsignal+0x110)[0x7f2946c1f0d0] 
>> [linux-lb8c:15830] [ 2] /lib64/libc.so.6(abort+0x151)[0x7f2946c206b1] 
>> [linux-lb8c:15830] [ 3] 
>> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(+0x6b9e5d1)[0x7f295b49e5d1] 
>> [linux-lb8c:15830] [ 4] 
>> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals5abortERKNS_13ExceptionBaseE+0x1a)[0x7f295b49edaf]
>>  
>> [linux-lb8c:15830] [ 5] 
>> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals11issue_errorINS_18StandardExceptions17ExcNotImplementedEEEvNS1_17ExceptionHandlingEPKciS7_S7_S7_T_+0x98)[0x7f2957373ea1]
>>  
>> [linux-lb8c:15830] [ 6] 
>> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(+0x3f38e23)[0x7f2958838e23] 
>> [linux-lb8c:15830] [ 7] 
>> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii11VectorTools7projectILi2ENS_16TrilinosWrappers3MPI6VectorELi2EEEvRKNS_10DoFHandlerIXT_EXT1_EEERKNS_16ConstraintMatrixERKNS_10QuadratureIXT_EEERKNS_8FunctionIXT1_ENT0_10value_typeEEERSH_bRKNSC_IX
>> miT_Li1b+0x2f)[0x7f295894906e] 
>> [linux-lb8c:15830] [ 8] 
>> ./main(_ZN6Step1521MinimalSurfaceProblemILi2EE3runEv+0xc08)[0x420d08] 
>> [linux-lb8c:15830] [ 9] ./main(main+0x3c)[0x414ad0] 
>> [linux-lb8c:15830] [10] 
>> /lib64/libc.so.6(__libc_start_main+0xea)[0x7f2946c09f4a] 
>> [linux-lb8c:15830] [11] ./main(_start+0x2a)[0x41477a] 
>> [linux-lb8c:15830] *** End of error message *** 
>> Abgebrochen (Speicherabzug geschrieben)
>>
>> when running in debug mode. It runs fine in release mode. Why does that 
>> happen for more than two components, and how can I fix/circumvent that? Or 
>> did I (again) forget something? 
>>
>> My minimal example is attached, the behaviour happens when setting 
>> NUM_COMPONENTS via 
>>
>> #define NUM_COMPONENTS 100
>>
>> to a value larger than 2.
>>
>>
>> Thank you!
>>
>>

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

Re: [deal.II] Solving time dependent heat equation with MPI (PETsc)

[Wolfgang Bangerth]

On 10/16/2017 08:55 AM, Mark Ma wrote:


 > So when you visualize the solution, the error is at the boundary
but it looks
 > correct in the interior?
 >
 > Yes, it is. Thereafter, the error at boundaries does propagate and then
 > interfere with the interior.

  From the pictures you posted in a follow-up, it looks like the boundary
values are actually correct, but that the problem starts one layer of cells
into the domain. Do I see this correctly?


Yes, exactly. And if I use the BC values as 0, then it seems OK now (see figs 
below). What if I want to apply some value other than 0, this problem seems 
annoying.


Yes, bugs are highly annoying indeed :-)

It almost looks to me like you're applying both Dirichlet values (via the 
ConstraintMatrix) and Neumann boundary values (via a boundary integral). But 
then I haven't taken a look at the code, so I can't really say for sure.


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.

Re: [deal.II] Solving time dependent heat equation with MPI (PETsc)

[Mark Ma]

在 2017年10月16日星期一 UTC+2下午4:31:25，Wolfgang Bangerth写道：
>
>
> > So when you visualize the solution, the error is at the boundary but 
> it looks 
> > correct in the interior? 
> > 
> > Yes, it is. Thereafter, the error at boundaries does propagate and then 
> > interfere with the interior. 
>
>  From the pictures you posted in a follow-up, it looks like the boundary 
> values are actually correct, but that the problem starts one layer of 
> cells 
> into the domain. Do I see this correctly? 
>
>
Yes, exactly. And if I use the BC values as 0, then it seems OK now (see 
figs below). What if I want to apply some value other than 0, this problem 
seems annoying.
 

>
> >  > Is there a simple way to update the RHS of old value using 
> something 
> > simple 
> >  > like vmult? 
> > 
> > What is the equation you are trying to implement? 
> > 
> > I am trying just the Heat equation, 
> > rho * C * dT/dt = nabla *(kappba nalbaT) + f(x,y,t); 
> > Boundary condition is, 
> > T(x,y,z,t) = 300K, at boundaries. 
> I meant the "update the RHS of old value" equation. 
>   
>

sorry, I mistook your idea. The RHS function is,
f(x,y,z,t) = I0 * exp(-((x-v*t)^2 + (y-v*t)^2 +(z-v*t)^2)/2c^2), which 
stands for absorption of the laser.





Boundary value = 0, initial values = 0


-- 
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] Solving time dependent heat equation with MPI (PETsc)

[Wolfgang Bangerth]

So when you visualize the solution, the error is at the boundary but it 
looks
correct in the interior?

Yes, it is. Thereafter, the error at boundaries does propagate and then 
interfere with the interior.


From the pictures you posted in a follow-up, it looks like the boundary 
values are actually correct, but that the problem starts one layer of cells 
into the domain. Do I see this correctly?




 > Is there a simple way to update the RHS of old value using something
simple
 > like vmult?

What is the equation you are trying to implement?

I am trying just the Heat equation,
rho * C * dT/dt = nabla *(kappba nalbaT) + f(x,y,t);
Boundary condition is,
T(x,y,z,t) = 300K, at boundaries.

I meant the "update the RHS of old value" equation.

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.

[deal.II] Re: Solving time dependent heat equation with MPI (PETsc)

[Mark Ma]

Projection of initial value, time 000








time 001





time 002










在 2017年10月10日星期二 UTC+2下午4:40:34，Mark Ma写道：
>
>
> Dear experienced deal.II users and developers, 
>
> I want to solve a heat equation in the time domain with distributed memory 
> using MPI, but the results are incorrect. In order to do so, I reference 
> tutorial step-23 for time updating method and step-40 for implementing MPI. 
> May I ask whether my boundary condition is right or not? Should we do 
> compress() after apply_boundary_values()? Thanks in advance!
>
> With best regards,
> Mark
>
>
>
>
>
> In PETsc or Trilinos, how could we set the boundary conditions and initial 
> conditions?
>
>
> Since old_solution_T would be read in other places (not shown here), it is 
> initialized with ghost cells; while old_solutuon_T_cal is only used for 
> written, it does not have ghost cells. see code as following,
> //
> old_solution_T_cal.reinit (locally_owned_dofs,mpi_communicator);
> old_solution_T.reinit 
> (locally_owned_dofs,locally_relevant_dofs,mpi_communicator);
> //
> ..
>
>
>   template 
>   void ThermalDiffusion::run ()
>   {
>
>
> setup_system();
> assemble_system();
>
> 
> VectorTools::project (dof_handler, constraints, QGauss(degree),
>   InitialValues_T(),
>   old_solution_T_cal);
>
> old_solution_T = old_solution_T_cal;
>
> LA::MPI::Vector tmp (locally_owned_dofs,mpi_communicator);
> LA::MPI::Vector forcing_terms (locally_owned_dofs,mpi_communicator);
>
> for (timestep_number=1, time=time_step;
>  time<= global_simulation_end_time;
>  time+=time_step, ++timestep_number)
>   {
> pcout << "Time step " << timestep_number
>   << " at t=" << time
>   << std::endl;
>
>
> //---
> //run to solve T
> //---
>
>   //
>   //time dependent
>
> //assign right hand side
> mass_matrix_T.vmult (system_rhs, old_solution_T_cal);
>
> laplace_matrix_T.vmult (tmp, old_solution_T_cal);
> system_rhs.add (-time_step * (1-theta), tmp);
>
> assemble_rhs_T (time);
> forcing_terms = dynamic_rhs_T;
> forcing_terms *= time_step * theta;
>
> assemble_rhs_T (time - time_step); 
> tmp = dynamic_rhs_T;
> forcing_terms.add (time_step*(1-theta),tmp);
> system_rhs.add (1,forcing_terms);
>
> //assign system matrix
> system_matrix.copy_from (mass_matrix_T);
> system_matrix.add (time_step * theta, laplace_matrix_T);
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> * {  BoundaryValues_Temperature 
> boundary_values_function_T;  //boundary_values_function.set_time 
> (time);  std::map 
> boundary_values_T;  VectorTools::interpolate_boundary_values 
> (dof_handler,
> BOUNDARY_NUM,
> boundary_values_function_T,   
>  
> boundary_values_T);  MatrixTools::apply_boundary_values 
> (boundary_values_T,  
> system_matrix,  
> solution_T,  
> system_rhs,  

Re: [deal.II] Solving time dependent heat equation with MPI (PETsc)

[Mark Ma]

> On 10/15/2017 12:41 PM, Mark Ma wrote: 
> > 
> > Now the projection of initial values (rewrite the code by manually 
> assemble 
> > the matrix and system_rhs and calculate) run OK, but the time updating 
> of T is 
> > not correct, same phenomenon appears. I believe this may arise from the 
> fact 
> > that direct using matrix vmult (i.e. 
> > | 
> > mass_matrix_T.vmult (system_rhs,old_solution_T_cal); 
> > | 
> > 
> > ) instead of assembling and distribute_local_to_global again may ignore 
> > eliminating the constraint points in matrix or vector, when using 
> > constriantMatrix and interporate_boundary_values to apply the boundary 
> > condition, I am now checking of this. 
>
> So when you visualize the solution, the error is at the boundary but it 
> looks 
> correct in the interior? 
>
> Yes, it is. Thereafter, the error at boundaries does propagate and then 
interfere with the interior. 

>
> > Is there a simple way to update the RHS of old value using something 
> simple 
> > like vmult? 
>
> What is the equation you are trying to implement? 
>
> I am trying just the Heat equation,
rho * C * dT/dt = nabla *(kappba nalbaT) + f(x,y,t);
Boundary condition is,
T(x,y,z,t) = 300K, at boundaries.
 

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

[deal.II] Re: Error when applying initial values to MPI::Vector in multiple dimensions

[Denis Davydov]

Or you may want to use interpolate if this is enough
https://www.dealii.org/developer/doxygen/deal.II/namespaceVectorTools.html#a05db6c8cebf924b417dd92f525efe3db

Regards,
Denis

On Friday, October 13, 2017 at 11:05:01 AM UTC+2, Maxi Miller wrote:
>
> I try to apply initial values to a vector defined as 
> *LinearAlgebraTrilinos::MPI::Vector 
> *using 
> VectorTools::project (dof_handler, hanging_node_constraints,
>  QGauss(fe.degree+1),
>  InitialValues(),
>  local_solution);
>
>
>
> When initializing the variable fe (as FESystem) with one or two 
> components, it works fine. For more than two components I get the error
>
>  
> An error occurred in line <1366> of file 
> <~/Downloads/dealii/include/deal.II/numerics/vector_tools.templates.h> in 
> function 
> void dealii::VectorTools::{anonymous}::project(const 
> dealii::Mapping&, const dealii::DoFHandler&, const 
> dealii::ConstraintMatrix&, const dealii::Quadrature&, const 
> dealii::Function&, VectorType&, bool, 
> const dealii
> ::Quadrature<(dim - 1)>&, bool) [with VectorType = 
> dealii::TrilinosWrappers::MPI::Vector; int dim = 2; typename 
> VectorType::value_type = double] 
> The violated condition was:  
> (dynamic_cast* > 
> (&(dof.get_triangulation()))==nullptr) 
> Additional information:  
> You are trying to use functionality in deal.II that is currently not 
> implemented. In many cases, this indicates that there simply didn't appear 
> much of a need for it, or that the author of the original code did not have 
> the time to implement a particular case. If you
>  hit this exception, it is therefore worth the time to look into the code to 
> find out whether you may be able to implement the missing functionality. If 
> you do, please consider providing a patch to the deal.II development sources 
> (see the deal.II website on how to contri
> bute). 
>  
> Stacktrace: 
> --- 
> #0  /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre:  
> #1  /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre: void 
> dealii::VectorTools::project<2, dealii::TrilinosWrappers::MPI::Vector, 
> 2>(dealii::DoFHandler<2, 2> const&, dealii::ConstraintMatrix const&, 
> dealii::Quadrature<2> const&, dealii::Function<2, 
> dealii::TrilinosWrappers::MPI
> ::Vector::value_type> const&, dealii::TrilinosWrappers::MPI::Vector&, bool, 
> dealii::Quadrature<(2)-(1)> const&, bool) 
> #2  ./main: Step15::MinimalSurfaceProblem<2>::run() 
> #3  ./main: main 
>  
>  
> [linux-lb8c:15830] *** Process received signal *** 
> [linux-lb8c:15830] Signal: Aborted (6) 
> [linux-lb8c:15830] Signal code:  (-6) 
> [linux-lb8c:15830] [ 0] /lib64/libpthread.so.0(+0x12270)[0x7f294a477270] 
> [linux-lb8c:15830] [ 1] /lib64/libc.so.6(gsignal+0x110)[0x7f2946c1f0d0] 
> [linux-lb8c:15830] [ 2] /lib64/libc.so.6(abort+0x151)[0x7f2946c206b1] 
> [linux-lb8c:15830] [ 3] 
> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(+0x6b9e5d1)[0x7f295b49e5d1] 
> [linux-lb8c:15830] [ 4] 
> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals5abortERKNS_13ExceptionBaseE+0x1a)[0x7f295b49edaf]
>  
> [linux-lb8c:15830] [ 5] 
> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals11issue_errorINS_18StandardExceptions17ExcNotImplementedEEEvNS1_17ExceptionHandlingEPKciS7_S7_S7_T_+0x98)[0x7f2957373ea1]
>  
> [linux-lb8c:15830] [ 6] 
> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(+0x3f38e23)[0x7f2958838e23] 
> [linux-lb8c:15830] [ 7] 
> /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii11VectorTools7projectILi2ENS_16TrilinosWrappers3MPI6VectorELi2EEEvRKNS_10DoFHandlerIXT_EXT1_EEERKNS_16ConstraintMatrixERKNS_10QuadratureIXT_EEERKNS_8FunctionIXT1_ENT0_10value_typeEEERSH_bRKNSC_IX
> miT_Li1b+0x2f)[0x7f295894906e] 
> [linux-lb8c:15830] [ 8] 
> ./main(_ZN6Step1521MinimalSurfaceProblemILi2EE3runEv+0xc08)[0x420d08] 
> [linux-lb8c:15830] [ 9] ./main(main+0x3c)[0x414ad0] 
> [linux-lb8c:15830] [10] 
> /lib64/libc.so.6(__libc_start_main+0xea)[0x7f2946c09f4a] 
> [linux-lb8c:15830] [11] ./main(_start+0x2a)[0x41477a] 
> [linux-lb8c:15830] *** End of error message *** 
> Abgebrochen (Speicherabzug geschrieben)
>
> when running in debug mode. It runs fine in release mode. Why does that 
> happen for more than two components, and how can I fix/circumvent that? Or 
> did I (again) forget something? 
>
> My minimal example is attached, the behaviour happens when setting 
> NUM_COMPONENTS via 
>
> #define NUM_COMPONENTS 100
>
> to a value larger than 2.
>
>
> 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.