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.
