Hi Jaekwang,
The error message generated by PETSc states that you are trying to insert
an invalid floating point number into the system matrix. This probably
means that you are using uninitialized values somewhere in the computation:
usually this mistake is caught by deal.II. Did you compile your deal.II
application in debug mode? If not, I suspect that running the program again
in debug mode will show exactly where something was not initialized
(deal.II checks this in many places).
Thanks,
David Wells
On Thursday, September 21, 2017 at 2:16:26 PM UTC-4, Jaekwang Kim wrote:
>
> Hi, all
>
> I was practicing using Petsc library through deal.ii wrappers,
>
> While, step 17 gives me clear explanation, I am having trouble in
> distributing hanging node constraints on Petsc system matrix,
>
> I was modifying step-9 hyperbolic equation solvers with PETSc matrix, but
> when I compile my code I run into error message like
>
> *[0]PETSC ERROR: --------------------- Error Message
> ------------------------------------*
>
> *[0]PETSC ERROR: Floating point exception!*
>
> *[0]PETSC ERROR: Inserting nan at matrix entry (0,0)!*
>
> *[0]PETSC ERROR:
> ------------------------------------------------------------------------*
>
> *[0]PETSC ERROR: Petsc Release Version 3.3.0, Patch 6, Mon Feb 11 12:26:34
> CST 2013 *
>
> *[0]PETSC ERROR: See docs/changes/index.html for recent updates.*
>
> *[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.*
>
> *[0]PETSC ERROR: See docs/index.html for manual pages.*
>
> *[0]PETSC ERROR:
> ------------------------------------------------------------------------*
>
> *[0]PETSC ERROR: ./advection on a dbg named golubh2 by jk12 Thu Sep 21
> 13:05:46 2017*
>
> *[0]PETSC ERROR: Libraries linked from
> /usr/local/apps/petsc/3.3-p6/intel/mvapich2-1.6-intel/dbg/lib*
>
> *[0]PETSC ERROR: Configure run at Wed Mar 13 14:51:28 2013*
>
> *[0]PETSC ERROR: Configure options
> --prefix=/usr/local/apps/petsc/3.3-p6/intel/mvapich2-1.6-intel/dbg
> --with-cc=mpicc --with-fc=mpif90 --with-cxx=mpicxx
> --with-shared-libraries=0 --with-mpi=1 --known-mpi-shared=1
> --with-blas-lapack-lib="[/usr/local/intel-11.1/mkl/lib/em64t/libmkl_intel_lp64.a,/usr/local/intel-11.1/mkl/lib/em64t/libmkl_sequential.a,/usr/local/intel-11.1/mkl/lib/em64t/libmkl_core.a]"*
>
> *[0]PETSC ERROR:
> ------------------------------------------------------------------------*
>
> *[0]PETSC ERROR: MatSetValues() line 1013 in
> /usr/local/apps/petsc/build/petsc-3.3-p6/src/mat/interface/matrix.c*
>
> *ERROR: Uncaught exception in MPI_InitFinalize on proc 0. Skipping
> MPI_Finalize() to avoid a deadlock.*
>
> I have checked where I got stuck in my code and found that they are
> generated when I distribute constraints to (yellow highlighted lines).
>
> I attached how I made constraints matrix at setup system.
>
> I would appreciate if any one can suggests where I got stuck...
>
> Thanks
>
> Jaekwang
>
>
>
>
>
> template <int dim>
>
> void AdvectionProblem<dim>::assemble_system_petsc ()
>
> {
>
> pcout << "assemble system A" << std::endl;
>
>
>
> const MappingQ<dim> mapping (degree);
>
> QGauss<dim> quadrature_formula(3);
>
> QGauss<dim-1> face_quadrature_formula(3);
>
>
>
> FEValues<dim> fe_values (mapping, fe, quadrature_formula,
>
> update_values | update_gradients |
>
> update_quadrature_points |
> update_JxW_values);
>
>
>
> FEFaceValues<dim> fe_face_values (mapping, fe,
> face_quadrature_formula,
>
> update_values |
> update_normal_vectors | update_gradients |
>
> update_quadrature_points |
> update_JxW_values);
>
>
>
> const AdvectionField<dim> advection_field; //Defined by Tensor
> Function
>
> const RightHandSide<dim> right_hand_side; //Defined by normal
> Vector Function
>
> const BoundaryValues<dim> boundary_values; //Defined by normal
> Vector Function
>
>
>
> const unsigned int dofs_per_cell = fe.dofs_per_cell;
>
> const unsigned int n_q_points = quadrature_formula.size();
>
> const unsigned int n_face_q_points =
> fe_face_values.get_quadrature().size();
>
>
>
> FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
>
> Vector<double> cell_rhs (dofs_per_cell);
>
>
>
> std::vector<types::global_dof_index> local_dof_indices
> (dofs_per_cell);
>
>
>
> std::vector<double> rhs_values (n_q_points);
>
> std::vector<Tensor<1,dim> > advection_directions (n_q_points);
>
> std::vector<double> face_boundary_values
> (n_face_q_points);
>
> std::vector<Tensor<1,dim> > face_advection_directions
> (n_face_q_points);
>
>
>
> double tau ;
>
>
>
> pcout << "assemble system B" << std::endl;
>
>
>
> typename DoFHandler<dim>::active_cell_iterator
>
> cell = dof_handler.begin_active(),
>
> endc = dof_handler.end();
>
> for (; cell!=endc; ++cell)
>
> if (cell->subdomain_id() == this_mpi_process)
>
> {
>
>
>
> pcout << "assemble system C" << std::endl;
>
>
>
> cell_matrix = 0;
>
> cell_rhs = 0;
>
>
>
> tau = 0.5 * cell->diameter () * cell->diameter ();
>
>
>
> advection_field.value_list
> (fe_values.get_quadrature_points(),
>
> advection_directions);
>
> right_hand_side.value_list
> (fe_values.get_quadrature_points(),
>
> rhs_values);
>
>
>
>
>
> // pcout << "assemble system D " << std::endl;
>
>
>
> for (unsigned int q_point=0; q_point<n_q_points;
> ++q_point)
>
> for (unsigned int i=0; i<dofs_per_cell; ++i)
>
> {
>
> for (unsigned int j=0; j<dofs_per_cell; ++j) //
> v-[i] , u [j]
>
>
>
> {
>
> //pcout << "assemble system E " << std::endl;
>
>
>
> cell_matrix(i,j) +=
>
> (
>
> // Advection Term
>
> fe_values.shape_value(i,q_point) * (
> advection_directions[q_point] * fe_values.shape_grad(j,q_point)) //(beta
> grad u)
>
>
>
> // Stabilization term (considering residual)
>
> + tau * ( advection_directions[q_point] *
> fe_values.shape_grad(i,q_point) )
>
> *(advection_directions[q_point] *
> fe_values.shape_grad(j,q_point)
>
>
>
> )
>
>
>
> ) * fe_values.JxW(q_point);
>
> }
>
>
>
> // pcout << "assemble system F " << std::endl;
>
>
>
> cell_rhs(i) += fe_values.shape_value(i,q_point) *
> rhs_values[q_point] * fe_values.JxW (q_point);
>
>
>
> }
>
>
>
> //Bonndary that has incoming flow...
>
>
>
>
>
> pcout << "assemble system G " << std::endl;
>
> for (unsigned int face=0;
> face<GeometryInfo<dim>::faces_per_cell; ++face)
>
> if (cell->face(face)->at_boundary())
>
> {
>
> pcout << "assemble system H " << std::endl;
>
>
>
> fe_face_values.reinit (cell, face);
>
>
>
> boundary_values.value_list
> (fe_face_values.get_quadrature_points(),
>
> face_boundary_values);
>
> advection_field.value_list
> (fe_face_values.get_quadrature_points(),
>
>
> face_advection_directions);
>
>
>
> for (unsigned int q_point=0;
> q_point<n_face_q_points; ++q_point)
>
> if (fe_face_values.normal_vector(q_point) *
>
> face_advection_directions[q_point]< 0)//this
> describes inflow case
>
> for (unsigned int i=0; i<dofs_per_cell;
> ++i)
>
> {
>
> for (unsigned int j=0;
> j<dofs_per_cell; ++j)
>
> cell_matrix(i,j) -=
> (face_advection_directions[q_point] *
>
>
> fe_face_values.normal_vector(q_point) *
>
>
> fe_face_values.shape_value(i,q_point) *
>
>
> fe_face_values.shape_value(j,q_point) *
>
>
> fe_face_values.JxW(q_point));
>
>
>
> cell_rhs(i) -=
> (face_advection_directions[q_point] *
>
>
> fe_face_values.normal_vector(q_point) *
>
>
> face_boundary_values[q_point] *
>
>
> fe_face_values.shape_value(i,q_point) *
>
>
> fe_face_values.JxW(q_point));
>
> }
>
> }// end of integrals on inflow boundary
>
>
>
>
> cell->get_dof_indices (local_dof_indices);
>
>
>
> pcout << "assemble system J " << std::endl;
>
>
>
>
>
>
> hanging_node_constraints.distribute_local_to_global(cell_matrix, cell_rhs,
>
> local_dof_indices,
>
> system_matrix_petsc,
> system_rhs_petsc);
>
>
>
>
>
> pcout << "assemble system K " << std::endl;
>
>
>
> }//end of cell iterator
>
>
>
>
>
> } // end of assemble_system_petsc
>
>
>
>
> template <int dim>
>
> void AdvectionProblem<dim>::setup_system_petsc ()
>
> {
>
> pcout <<" setup_system_petsc started" << std::endl;
>
>
>
> //Split triangulation into n-processes
>
> GridTools::partition_triangulation (n_mpi_processes,
> triangulation);
>
> dof_handler.distribute_dofs (fe);
>
>
>
> // Now you have multiple sub-domain on each proceses,
>
> // You need to renumbering Dof according to this data structure
>
> DoFRenumbering::subdomain_wise (dof_handler);
>
>
>
> hanging_node_constraints.clear ();
>
> DoFTools::make_hanging_node_constraints (dof_handler,
>
> hanging_node_constraints);
>
> hanging_node_constraints.close ();
>
>
>
>
>
> // therefore, kill them
>
> DynamicSparsityPattern dsp(dof_handler.n_dofs(),
> dof_handler.n_dofs());
>
> DoFTools::make_sparsity_pattern(dof_handler,
>
> dsp,
>
> hanging_node_constraints,
>
> /*keep_constrained_dofs = */ true
> );
>
>
>
> sparsity_pattern.copy_from (dsp);
>
>
>
> // prepare matrix and vectors that can be used with MPI
>
> const types::global_dof_index n_local_dofs =
> DoFTools::count_dofs_with_subdomain_association (dof_handler,
> this_mpi_process);
>
>
>
>
>
> system_matrix_petsc.reinit (mpi_communicator,
>
> dof_handler.n_dofs(),
>
> dof_handler.n_dofs(),
>
> n_local_dofs,
>
> n_local_dofs,
>
> dof_handler.max_couplings_between_dofs());
>
>
>
> solution_petsc.reinit (mpi_communicator, dof_handler.n_dofs(),
> n_local_dofs);
>
> system_rhs_petsc.reinit (mpi_communicator, dof_handler.n_dofs(),
> n_local_dofs);
>
>
>
>
>
> //generate exact solution
>
> exact_solution.reinit (dof_handler.n_dofs());
>
>
> VectorTools::interpolate(dof_handler,Exact_Solution<dim>(),exact_solution);
>
>
>
> }
>
>
--
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 [email protected].
For more options, visit https://groups.google.com/d/optout.
