Re: [deal.II] Re: Porting tutorials to PETSc from Trilinos
Thanks Daniel, now it runs. Of course it won't converge, lacking preconditiones, but this is for another question. Thank you! Franco -- 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/f26d4349-7889-41d6-8349-71a33fa31f71%40googlegroups.com.
Re: [deal.II] Re: Porting tutorials to PETSc from Trilinos
Franco, If you forget about the static_cast in the error message just above >> you can actually see what happens: >> The violated condition was: >> global_matrix.m() == global_vector.size() >> >> > > I am lost in this, I have made some stupid error here but I don't see how. > I am posting the code here: > > https://gist.github.com/fmilicchio/2dac430d7c071f98b4f66d79d2101b04 > > As far as I can tell, I've initialized every sparsity pattern, matrix, and > vector as they were with Trilinos, but I am surely seeing the code I > intended to write, not the code I have actually written. > You are not initializing the stokes vectors appropriately. Running in a debugger should have told you that there was a problem with the size of stokes_rhs. In fact, stokes_solution.reinit(n_u+n_p) creates a vector with n_u+n_p empty blocks. You want to do something like std::vector block_sizes(2); block_sizes[0] = n_u; block_sizes[1] = n_p; stokes_solution.reinit(block_sizes); instead. Then, you also need to have a look how to initialize the sparsity pattern for the stokes_preconditioner correctly (if you want to use that matrix at all). Best, Daniel -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/de31075a-c077-4ca5-b95e-84c0b5dfebd3%40googlegroups.com.
Re: [deal.II] Re: Porting tutorials to PETSc from Trilinos
> No unfortunately we don't have such things. However, most of the new > effort on linear algebra is in the LinearAlgebra namespace. > > Thanks, Bruno, I've started looking into that namespace. There are several options to make it work in a coherent way, I'll ask for comments on these solutions as I figured all classes out (and how to use MPI seamlessly, at first I am thinking about tag dispatching to avoid breaking legacy codes). > > I don't know why this error is so unreadable but if you look at line > 3781 in include/deal.II/lac/affine_constraints.templates.h you will > see that the assert checks that the global matrix has as many rows has > the global vector. To me it looks like the problem is that the global > vector as a size of zero. > > If you forget about the static_cast in the error message just above > you can actually see what happens: > The violated condition was: > global_matrix.m() == global_vector.size() > > I am lost in this, I have made some stupid error here but I don't see how. I am posting the code here: https://gist.github.com/fmilicchio/2dac430d7c071f98b4f66d79d2101b04 As far as I can tell, I've initialized every sparsity pattern, matrix, and vector as they were with Trilinos, but I am surely seeing the code I intended to write, not the code I have actually written. Do you see any stupid error at first sight? Thank you! Franco -- 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/ce549f8c-9a44-462f-a66f-fcfab2794b3e%40googlegroups.com.
Re: [deal.II] Re: Porting tutorials to PETSc from Trilinos
Franco Le ven. 19 juil. 2019 à 10:04, Franco Milicchio a écrit : > Is there a small documentation that I can read? For instance, I see there are > several namespaces like LinearAlgebra, PETSc/Epetra/Trilinos/CUDAWrappers, > LinearAlgebraDealii and so on, that might be restructured without > compromising old codes, adding some uniformity to classes. No unfortunately we don't have such things. However, most of the new effort on linear algebra is in the LinearAlgebra namespace. > > An error occurred in line <3781> of file > > in function > void dealii::AffineConstraints::distribute_local_to_global(const > FullMatrix &, const Vector &, const std::vector &, > MatrixType &, VectorType &, bool, std::integral_constant) const > [number = double, MatrixType = dealii::BlockSparseMatrix, VectorType > = dealii::BlockVector] > The violated condition was: > static_cast std::common_type decltype(global_vector.size())>::type)>::type>(global_matrix.m()) == > static_cast std::common_type decltype(global_vector.size())>::type)>::type>(global_vector.size()) > Additional information: I don't know why this error is so unreadable but if you look at line 3781 in include/deal.II/lac/affine_constraints.templates.h you will see that the assert checks that the global matrix has as many rows has the global vector. To me it looks like the problem is that the global vector as a size of zero. If you forget about the static_cast in the error message just above you can actually see what happens: The violated condition was: global_matrix.m() == global_vector.size() 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. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/CAGVt9eM%2Bg5HhuJOyQ3%3D-D__EbUyf7tZS5ZqyTu32rxVN1S-vXw%40mail.gmail.com.
Re: [deal.II] Re: Porting tutorials to PETSc from Trilinos
> > Thanks for your answer. By the way, would a wrapper for the Eigen > library be appealing to anyone? > I am sure it would be useful! > > Good to know, Bruno. I will start looking into those wrappers and I hope some student of mine will be helping (fingers crossed). Is there a small documentation that I can read? For instance, I see there are several namespaces like LinearAlgebra, PETSc/Epetra/Trilinos/CUDAWrappers, LinearAlgebraDealii and so on, that might be restructured without compromising old codes, adding some uniformity to classes. > > It looks like the sparsity pattern does not exist anymore. The matrix > object does not keep a copy of the sparsity pattern, so the sparsity > pattern needs to exist as long as the matrix that uses it is alive. > > Thanks, that did the trick! The problem now moved to the AffineConstraints (which I didn't touch) complaining about dimensions: An error occurred in line <3781> of file in function void dealii::AffineConstraints::distribute_local_to_global(const FullMatrix &, const Vector &, const std::vector &, MatrixType &, VectorType &, bool, std::integral_constant) const [number = double, MatrixType = dealii::BlockSparseMatrix, VectorType = dealii::BlockVector] The violated condition was: static_cast::type)>::type>(global_matrix.m()) == static_cast::type)>::type>(global_vector.size()) Additional information: Dimension 2467 not equal to 0. Stacktrace: --- #0 2 libdeal_II.g.9.1.0.dylib0x0001164eb3a2 _ZNK6dealii17AffineConstraintsIdE26distribute_local_to_globalINS_17BlockSparseMatrixIdEENS_11BlockVectorIdvRKNS_10FullMatrixIdEERKNS_6VectorIdEERKNSt3__16vectorIjNSF_9allocatorIjRT_RT0_bNSF_17integral_constantIbLb1EEE + 1570: 2 libdeal_II.g.9.1.0.dylib0x0001164eb3a2 _ZNK6dealii17AffineConstraintsIdE26distribute_local_to_globalINS_17BlockSparseMatrixIdEENS_11BlockVectorIdvRKNS_10FullMatrixIdEERKNS_6VectorIdEERKNSt3__16vectorIjNSF_9allocatorIjRT_RT0_bNSF_17integral_constantIbLb1EEE #1 3 step-31 0x000108d5b641 _ZN6Step3121BoussinesqFlowProblemILi2EE22assemble_stokes_systemEv + 1953: 3 step-31 0x000108d5b641 _ZN6Step3121BoussinesqFlowProblemILi2EE22assemble_stokes_systemEv #2 4 step-31 0x000108d4e468 _ZN6Step3121BoussinesqFlowProblemILi2EE3runEv + 504: 4 step-31 0x000108d4e468 _ZN6Step3121BoussinesqFlowProblemILi2EE3runEv #3 5 step-31 0x000108d4de6d main + 125: 5 step-31 0x000108d4de6d main #4 6 libdyld.dylib 0x7fff6bbd93d5 start + 1: 6 libdyld.dylib 0x7fff6bbd93d5 start #5 7 ??? 0x0001 0x0 + 1: 7 ??? 0x0001 0x0 Signal: SIGABRT (signal SIGABRT) on this line (initial timestep): if (rebuild_stokes_matrix == true) stokes_constraints.distribute_local_to_global(local_matrix, local_rhs, local_dof_indices, stokes_matrix, stokes_rhs); The stokes_constraints.distribute_local_to_global apparently doesn't like how I implemented the sparsity patterns to build the matrices. All I did here is to add two sparsity patterns, one for the Stokes matrix, and one for the temperature: stokes_matrix_sp.reinit(2, 2); // probably unnecessary stokes_matrix_sp.copy_from(dsp); stokes_matrix.reinit(stokes_matrix_sp); //stokes_matrix.reinit(dsp); BlockSparsityPattern stokes_matrix_sp; // the definition I am probably doing something terribly stupid here... Thanks for any help! Franco -- 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/a7e083cb-8448-4840-aa0d-acb73e1207df%40googlegroups.com.
Re: [deal.II] Re: Porting tutorials to PETSc from Trilinos
> > Thanks for your answer. By the way, would a wrapper for the Eigen > library be appealing to anyone? > I am sure it would be useful! > > Good to know, Bruno. I will start looking into those wrappers and I hope some student of mine will be helping (fingers crossed). Is there a small documentation that I can read? For instance, I see there are several namespaces like LinearAlgebra, PETSc/Epetra/Trilinos/CUDAWrappers, LinearAlgebraDealii and so on, that might be restructured without compromising old codes, adding some uniformity to classes. > > It looks like the sparsity pattern does not exist anymore. The matrix > object does not keep a copy of the sparsity pattern, so the sparsity > pattern needs to exist as long as the matrix that uses it is alive. > > Thanks, that did the trick! The problem now moved to the AffineConstraints (which I didn't touch) complaining about dimensions: An error occurred in line <3781> of file in function void dealii::AffineConstraints::distribute_local_to_global(const FullMatrix &, const Vector &, const std::vector &, MatrixType &, VectorType &, bool, std::integral_constant) const [number = double, MatrixType = dealii::BlockSparseMatrix, VectorType = dealii::BlockVector] The violated condition was: static_cast::type)>::type>(global_matrix.m()) == static_cast::type)>::type>(global_vector.size()) Additional information: Dimension 2467 not equal to 0. Stacktrace: --- #0 2 libdeal_II.g.9.1.0.dylib0x0001164eb3a2 _ZNK6dealii17AffineConstraintsIdE26distribute_local_to_globalINS_17BlockSparseMatrixIdEENS_11BlockVectorIdvRKNS_10FullMatrixIdEERKNS_6VectorIdEERKNSt3__16vectorIjNSF_9allocatorIjRT_RT0_bNSF_17integral_constantIbLb1EEE + 1570: 2 libdeal_II.g.9.1.0.dylib0x0001164eb3a2 _ZNK6dealii17AffineConstraintsIdE26distribute_local_to_globalINS_17BlockSparseMatrixIdEENS_11BlockVectorIdvRKNS_10FullMatrixIdEERKNS_6VectorIdEERKNSt3__16vectorIjNSF_9allocatorIjRT_RT0_bNSF_17integral_constantIbLb1EEE #1 3 step-31 0x000108d5b641 _ZN6Step3121BoussinesqFlowProblemILi2EE22assemble_stokes_systemEv + 1953: 3 step-31 0x000108d5b641 _ZN6Step3121BoussinesqFlowProblemILi2EE22assemble_stokes_systemEv #2 4 step-31 0x000108d4e468 _ZN6Step3121BoussinesqFlowProblemILi2EE3runEv + 504: 4 step-31 0x000108d4e468 _ZN6Step3121BoussinesqFlowProblemILi2EE3runEv #3 5 step-31 0x000108d4de6d main + 125: 5 step-31 0x000108d4de6d main #4 6 libdyld.dylib 0x7fff6bbd93d5 start + 1: 6 libdyld.dylib 0x7fff6bbd93d5 start #5 7 ??? 0x0001 0x0 + 1: 7 ??? 0x0001 0x0 Signal: SIGABRT (signal SIGABRT) on this line (initial timestep): if (rebuild_stokes_matrix == true) stokes_constraints.distribute_local_to_global(local_matrix, local_rhs, local_dof_indices, stokes_matrix, stokes_rhs); The stokes_constraints.distribute_local_to_global apparently doesn't like how I implemented the sparsity patterns to build the matrices. All I did here is to add two sparsity patterns, one for the Stokes matrix, and one for the temperature: stokes_matrix_sp.reinit(2, 2); // probably unnecessary stokes_matrix_sp.copy_from(dsp); stokes_matrix.reinit(stokes_matrix_sp); //stokes_matrix.reinit(dsp); BlockSparsityPattern stokes_matrix_sp; // the definition stokes_matrix_sp.copy_from(dsp); stokes_matrix.reinit(stokes_matrix_sp); //stokes_matrix.reinit(dsp); I am probably doing something terribly stupid here... Thanks for any help! Franco -- 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/3bbde225-e547-43ab-abf3-2ef53dd09559%40googlegroups.com.
Re: [deal.II] Re: Porting tutorials to PETSc from Trilinos
Franco, Le jeu. 18 juil. 2019 à 10:27, Franco Milicchio a écrit : > Thanks for your answer. By the way, would a wrapper for the Eigen library be > appealing to anyone? I am sure it would be useful! > > An error occurred in line <382> of file > > in function > T *dealii::SmartPointer dealii::SparseMatrix >::operator->() const [T = const > dealii::SparsityPattern, P = dealii::SparseMatrix] > The violated condition was: > pointed_to_object_is_alive > Additional information: > The object pointed to is not valid anymore. It looks like the sparsity pattern does not exist anymore. The matrix object does not keep a copy of the sparsity pattern, so the sparsity pattern needs to exist as long as the matrix that uses it is alive. 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. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/CAGVt9eNGQNQt9qe389vk8ZWm2e02z_T8_cQ53_d%2B5KQ9%3DdX-jw%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] Re: Porting tutorials to PETSc from Trilinos
Thanks for your answer. By the way, would a wrapper for the Eigen library be appealing to anyone? And back to my lowly question :) I was able to port the Step 31 code to plain deal.II without PETSc, but I have an error at runtime. At the moment I removed the preconditioners and just used a PreconditionIdentity to make it compile and link, I will deal with the right preconditioner later. Apparently I have initialized something wrong, as the assignment to the matrix if (rebuild_stokes_matrix == true) stokes_matrix = 0; makes the program die. Number of active cells: 256 (on 5 levels) Number of degrees of freedom: 3556 (n_u 2178 + n_p 289 + n_T 1089) Timestep 0: t=0 Assembling...stokes_matrix 2467 x 2467 blocks 2 x 2 An error occurred in line <382> of file in function T *dealii::SmartPointer >::operator->() const [T = const dealii::SparsityPattern, P = dealii::SparseMatrix] The violated condition was: pointed_to_object_is_alive Additional information: The object pointed to is not valid anymore. Stacktrace: --- #0 2 libdeal_II.g.9.1.0.dylib0x0001174e003e _ZNK6dealii12SmartPointerIKNS_15SparsityPatternENS_12SparseMatrixIdEEEptEv + 206: 2 libdeal_II.g.9.1.0.dylib0x0001174e003e _ZNK6dealii12SmartPointerIKNS_15SparsityPatternENS_12SparseMatrixIdEEEptEv #1 3 libdeal_II.g.9.1.0.dylib0x00011abc65e5 _ZN6dealii12SparseMatrixIdEaSEd + 69: 3 libdeal_II.g.9.1.0.dylib 0x00011abc65e5 _ZN6dealii12SparseMatrixIdEaSEd #2 4 step-31 0x00010ada5f3a _ZN6dealii17BlockSparseMatrixIdEaSEd + 90: 4 step-31 0x00010ada5f3a _ZN6dealii17BlockSparseMatrixIdEaSEd #3 5 step-31 0x00010ada0e77 _ZN6Step3121BoussinesqFlowProblemILi2EE22assemble_stokes_systemEv + 247: 5 step-31 0x00010ada0e77 _ZN6Step3121BoussinesqFlowProblemILi2EE22assemble_stokes_systemEv #4 6 step-31 0x00010ad944b8 _ZN6Step3121BoussinesqFlowProblemILi2EE3runEv + 504: 6 step-31 0x00010ad944b8 _ZN6Step3121BoussinesqFlowProblemILi2EE3runEv #5 7 step-31 0x00010ad93ebd main + 125: 7 step-31 0x00010ad93ebd main #6 8 libdyld.dylib 0x7fff6bbd93d5 start + 1: 8 libdyld.dylib 0x7fff6bbd93d5 start I have checked the sizes as you can see, and they seem to be properly initialized. So my guess is that I am doing something wrong with sparsity patterns, but I don't know how: BlockSparsityPattern sp(2, 2); sp.copy_from(dsp); stokes_matrix.reinit(sp); //stokes_matrix.reinit(dsp); // ... SparsityPattern sp(n_T, n_T); sp.copy_from(dsp); temperature_matrix.reinit(sp); //temperature_matrix.reinit(dsp); temperature_mass_matrix.reinit(sp); //temperature_mass_matrix.reinit(temperature_matrix); temperature_stiffness_matrix.reinit(sp); //temperature_stiffness_matrix.reinit(temperature_matrix); The code I've modified is available in full here: https://gist.github.com/fmilicchio/2dac430d7c071f98b4f66d79d2101b04 Any hints? Thank you very much for your help! Franco -- 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/49a9517b-cb68-4dfd-96fb-0e1ed53ec1e7%40googlegroups.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] Re: Porting tutorials to PETSc from Trilinos
Franco, On Tuesday, July 16, 2019 at 3:54:34 AM UTC-4, Franco Milicchio wrote: > > > Can I ask why there are some discrepancies between those two wrappers? > Apart from an obvious one, that they're different libraries, but uniformity > may help a lot IMO. > Because the wrappers have been written over time by different people who were not aware of how exactly the other wrapper works (most people use either Trilinos or PETSc not both). We are trying to get the interfaces to be more similar but it is very difficult because any change in the interface will break many codes. 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. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/533ea4e2-44fd-43f0-8bd7-93b77bdbd6ed%40googlegroups.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] Re: Porting tutorials to PETSc from Trilinos
>
>
> step-40 tries to allow switching between Trilinos and PETSc as easy as
> possible by using the deal.II/lac/generic_linear_algebra.h header and the
> defining the LA::MPI namespace.
>
Thanks for your answer, Daniel. Even with that example I find it a little
harder than I expected to use classes from PETSc instead of Trilinos. For
instance:
const LinearSolvers::InverseMatrix
mp_inverse(stokes_preconditioner_matrix.block(1, 1),
*Mp_preconditioner);
Cannot compile, I am trying to use a
static_cast(stokes_preconditioner_matrix.block(1,
1)) now, but I don't think it will be a good idea.
Another weak point in my porting is
template
BlockSchurPreconditioner::
BlockSchurPreconditioner(
const PETScWrappers::MPI::BlockSparseMatrix ,
const InverseMatrix ,
const PreconditionerTypeA &Apreconditioner)
: stokes_matrix()
, m_inverse()
, a_preconditioner(Apreconditioner)
, tmp(complete_index_set(stokes_matrix->block(1, 1).m()))
{}
where tmp (a mutable PETScWrappers::MPI::Vector) cannot be initialized.
Adding MPI_COMM_WORLD make it compile, though.
Again with matrices, I have made the following changes:
//temperature_matrix.copy_from(temperature_mass_matrix);
temperature_matrix = 0;
temperature_matrix.add(1.0, temperature_mass_matrix);
I hope the change works...
Can I ask why there are some discrepancies between those two wrappers?
Apart from an obvious one, that they're different libraries, but uniformity
may help a lot IMO.
>
>> [...]
>> std::shared_ptr Amg_preconditioner;
>>
>
> This looks inappropriate. Probably it should rather be
> PETScWrappers::PreconditionBoomerAMG;
>
Yes, I know :) I am first trying to make it compile, then I will try to use
a good preconditioner (Boomer requires Hypre, and it needs MPI, so right
now it is not an option).
Thanks!
Franco
--
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/0bce36bd-19a1-4f62-8be8-75eb3965e7ce%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
[deal.II] Re: Porting tutorials to PETSc from Trilinos
Franco, > [...] > As I didn't find a 1:1 correspondence between the two namespaces, I've > resorted to remove Trilinos classes as best as I could, for instance, > TrilinosWrappers::BlockSparseMatrix exists but in PETSc there is only > PETScWrappers::MPI::BlockSparseMatrix that I know of. And the two seem to > be not exactly equivalent. > step-40 tries to allow switching between Trilinos and PETSc as easy as possible by using the deal.II/lac/generic_linear_algebra.h header and the defining the LA::MPI namespace. > [...] > std::shared_ptr Amg_preconditioner; > This looks inappropriate. Probably it should rather be PETScWrappers::PreconditionBoomerAMG; Best, Daniel -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/bdb2c6c2-c7d2-4ab6-bed6-9941740297df%40googlegroups.com. For more options, visit https://groups.google.com/d/optout.
