Hello all
I am trying to rewrite step-43 in parallel follwing paths of step-32. I
want just first to solve the darcy system then after I will move on solving
the saturation equation in time. While trying to slove the darcy system I
am having this errror message below:
The trace back is easy to follow. As there is no additional information, I
am just stuck.
Any hint dealii creator? :)
An error occurred in line <3635> of file
</home/franckm/deal.ii-candi/tmp/unpack/deal.II-v8.5.0/source/fe/fe_values.cc>
in function
void dealii::FEValues<dim, spacedim>::initialize(dealii::UpdateFlags)
[with int dim = 2; int spacedim = 2]
The violated condition was:
(update_flags & update_normal_vectors) == false
Additional information:
(none)
Stacktrace:
-----------
#0 /home/franckm/deal.ii-candi/deal.II-v8.5.0/lib/libdeal_II.g.so.8.5.0:
dealii::FEValues<2, 2>::initialize(dealii::UpdateFlags)
#1 /home/franckm/deal.ii-candi/deal.II-v8.5.0/lib/libdeal_II.g.so.8.5.0:
dealii::FEValues<2, 2>::FEValues(dealii::FiniteElement<2, 2> const&,
dealii::Quadrature<2> const&, dealii::UpdateFlags)
#2 ./darcyparl:
Step43::Assembly::Scratch::DarcyPreconditioner<2>::DarcyPreconditioner(dealii::FiniteElement<2,
2> const&, dealii::Quadrature<2> const&, dealii::UpdateFlags,
dealii::FiniteElement<2, 2> const&, dealii::UpdateFlags)
#3 ./darcyparl:
Step43::Assembly::Scratch::DarcySystem<2>::DarcySystem(dealii::FiniteElement<2,
2> const&, dealii::Quadrature<2> const&, dealii::Quadrature<1> const&,
dealii::UpdateFlags, dealii::FiniteElement<2, 2> const&,
dealii::UpdateFlags)
#4 ./darcyparl: Step43::TwoPhaseFlowProblem<2>::assemble_darcy_system()
#5 ./darcyparl: Step43::TwoPhaseFlowProblem<2>::solve()
#6 ./darcyparl: Step43::TwoPhaseFlowProblem<2>::run()
#7 ./darcyparl: main
and far down....
[franckm-Inspiron-3543:13997] *** Process received signal ***
[franckm-Inspiron-3543:13997] Signal: Aborted (6)
[franckm-Inspiron-3543:13997] Signal code: (-6)
[franckm-Inspiron-3543:13996] *** Process received signal ***
[franckm-Inspiron-3543:13996] Signal: Aborted (6)
[franckm-Inspiron-3543:13996] Signal code: (-6)
......
[franckm-Inspiron-3543:13998] *** Process received signal ***
[franckm-Inspiron-3543:13998] Signal: Aborted (6)
[franckm-Inspiron-3543:13998] Signal code: (-6)
[franckm-Inspiron-3543:13996] [ 0]
/lib/x86_64-linux-gnu/libpthread.so.0(+0x10330) [0x7f226587f330]
[franckm-Inspiron-3543:13996] [ 1]
/lib/x86_64-linux-gnu/libc.so.6(gsignal+0x37) [0x7f2264ed2c37]
[franckm-Inspiron-3543:13996] [ 2]
/lib/x86_64-linux-gnu/libc.so.6(abort+0x148) [0x7f2264ed6028]
[franckm-Inspiron-3543:13996] [ 3]
/home/franckm/deal.ii-candi/deal.II-v8.5.0/lib/libdeal_II.g.so.8.5.0(_ZN6dealii18deal_II_exceptions9internals5abortERKNS_13ExceptionBaseEb+0x3b)
[0x7f226cf4d4da]
[franckm-Inspiron-3543:13996] [ 4]
/home/franckm/deal.ii-candi/deal.II-v8.5.0/lib/libdeal_II.g.so.8.5.0(_ZN6dealii18deal_II_exceptions9internals11issue_errorINS_12FEValuesBaseILi2ELi2EE20ExcInvalidUpdateFlagEEEvNS1_17ExceptionHandlingEPKciS8_S8_S8_T_+0x30)
[0x7f226b275b40]
[franckm-Inspiron-3543:13996] [ 5]
/home/franckm/deal.ii-candi/deal.II-v8.5.0/lib/libdeal_II.g.so.8.5.0(_ZN6dealii8FEValuesILi2ELi2EE10initializeENS_11UpdateFlagsE+0x95)
[0x7f226b557ef5]
[franckm-Inspiron-3543:13996] [ 6]
/home/franckm/deal.ii-candi/deal.II-v8.5.0/lib/libdeal_II.g.so.8.5.0(_ZN6dealii8FEValuesILi2ELi2EEC1ERKNS_13FiniteElementILi2ELi2EEERKNS_10QuadratureILi2EEENS_11UpdateFlagsE+0x62)
[0x7f226b599ab2]
here my assembly and local_assemble_darcy_system ect.
namespace Assembly
{
namespace Scratch
{
template <int dim>
struct DarcyPreconditioner
{
DarcyPreconditioner (const FiniteElement<dim> &darcy_fe,
const Quadrature<dim> &darcy_quadrature,
const UpdateFlags darcy_update_flags,
const FiniteElement<dim> &saturation_fe,
const UpdateFlags
saturation_update_flags); // constructor
DarcyPreconditioner (const DarcyPreconditioner &data); // coppy
operation
FEValues<dim> darcy_fe_values;
FEValues<dim> saturation_fe_values;
std::vector<double> saturation_values;
std::vector<Tensor<2,dim>> k_inverse_values;
std::vector<Tensor<1,dim> > grad_phi_p;
std::vector<Tensor<1,dim> > phi_u;
};
template <int dim>
DarcyPreconditioner<dim>::
DarcyPreconditioner (const FiniteElement<dim> &darcy_fe,
const Quadrature<dim> &darcy_quadrature,
const UpdateFlags darcy_update_flags,
const FiniteElement<dim> &saturation_fe,
const UpdateFlags
saturation_update_flags)
:
darcy_fe_values (darcy_fe, darcy_quadrature, darcy_update_flags),
saturation_fe_values (saturation_fe, darcy_quadrature,
saturation_update_flags),
saturation_values (darcy_quadrature.size()),
k_inverse_values (darcy_quadrature.size()),
grad_phi_p (darcy_fe.dofs_per_cell),
phi_u (darcy_fe.dofs_per_cell)
{}
template <int dim>
DarcyPreconditioner<dim>::
DarcyPreconditioner (const DarcyPreconditioner &scratch)
:
darcy_fe_values (scratch.darcy_fe_values.get_fe(),
scratch.darcy_fe_values.get_quadrature(),
scratch.darcy_fe_values.get_update_flags()),
saturation_fe_values (scratch.saturation_fe_values.get_fe(),
scratch.saturation_fe_values.get_quadrature(),
scratch.saturation_fe_values.get_update_flags()),
saturation_values (scratch.saturation_values),
k_inverse_values (scratch.k_inverse_values),
grad_phi_p (scratch.grad_phi_p),
phi_u (scratch.phi_u)
{}
// The next one is the scratch object used for the assembly of the
full
// Stokes system. Observe that we derive the StokesSystem scratch
class
// from the StokesPreconditioner class above. We do this because all
the
// objects that are necessary for the assembly of the preconditioner
are
// also needed for the actual matrix system and right hand side, plus
// some extra data. This makes the program more compact. Note also
that
// the assembly of the Stokes system and the temperature right hand
side
// further down requires data from temperature and velocity,
// respectively, so we actually need two FEValues objects for those
two
// cases.
template <int dim>
struct DarcySystem : public DarcyPreconditioner<dim>
{
DarcySystem (const FiniteElement<dim> &darcy_fe,
const Quadrature<dim> &darcy_quadrature,
const Quadrature<dim-1> &darcy_face_quadrature,
const UpdateFlags darcy_update_flags,
const FiniteElement<dim> &saturation_fe,
const UpdateFlags saturation_update_flags);
DarcySystem (const DarcySystem<dim> &data);
FEFaceValues<dim> darcy_fe_face_values;
//FEValues<dim> saturation_fe_values;
std::vector<double> pressure_rhs_values;
std::vector<double> boundary_values ;
//std::vector<Tensor<2,dim> > k_inverse_values;
//std::vector<double> saturation_values;
std::vector<Tensor<1,dim> > phi_i_u;
std::vector<double> div_phi_u;
std::vector<double> phi_p;
};
template <int dim>
DarcySystem<dim>::
DarcySystem (const FiniteElement<dim> &darcy_fe,
const Quadrature<dim> &darcy_quadrature,
const Quadrature<dim-1> &darcy_face_quadrature,
const UpdateFlags darcy_update_flags,
const FiniteElement<dim> &saturation_fe,
const UpdateFlags saturation_update_flags)
:
DarcyPreconditioner<dim> (darcy_fe, darcy_quadrature,
darcy_update_flags,
saturation_fe,
saturation_update_flags),
darcy_fe_face_values (darcy_fe, darcy_face_quadrature,
darcy_update_flags),
//darcy_fe_face_values (darcy_fe, darcy_face_quadrature,
update_normal_vectors),
pressure_rhs_values (darcy_quadrature.size()),
boundary_values (darcy_face_quadrature.size()),
//saturation_values (darcy_quadrature.size()),
phi_i_u (darcy_fe.dofs_per_cell),
div_phi_u (darcy_fe.dofs_per_cell),
phi_p (darcy_fe.dofs_per_cell)
{}
template <int dim>
DarcySystem<dim>::
DarcySystem (const DarcySystem<dim> &scratch)
:
DarcyPreconditioner<dim> (scratch),
darcy_fe_face_values (scratch.darcy_fe_face_values.get_fe(),
scratch.darcy_fe_face_values.get_quadrature(),
scratch.darcy_fe_face_values.get_update_flags()),
phi_i_u (scratch.phi_i_u),
div_phi_u (scratch.div_phi_u),
phi_p (scratch.phi_p),
pressure_rhs_values (scratch.pressure_rhs_values),
boundary_values (scratch.boundary_values)
{}
template <int dim>
void
TwoPhaseFlowProblem<dim>::
local_assemble_darcy_preconditioner (const typename
DoFHandler<dim>::active_cell_iterator &cell,
Assembly::Scratch::DarcyPreconditioner<dim> &scratch,
Assembly::CopyData::DarcyPreconditioner<dim> &data)
{
//std::cout << " building darcy preconditioner..." << std::endl;
//darcy_preconditioner_matrix = 0;
/*const QGauss<dim> quadrature_formula(darcy_degree+2);
FEValues<dim> darcy_fe_values (darcy_fe, quadrature_formula,
update_JxW_values |
update_values |
update_gradients |
update_quadrature_points);*/
/*FEValues<dim> saturation_fe_values (saturation_fe,
quadrature_formula,
update_values);*/
//scratch.stokes_fe_values.get_fe().dofs_per_cell;
const unsigned int dofs_per_cell =
scratch.darcy_fe_values.get_fe().dofs_per_cell;
const unsigned int n_q_points =
scratch.darcy_fe_values.n_quadrature_points;
//std::vector<Tensor<2,dim> > k_inverse_values (n_q_points);
//std::vector<double> saturation_values (n_q_points);
//FullMatrix<double> local_matrix (dofs_per_cell,
dofs_per_cell);
//std::vector<types::global_dof_index> local_dof_indices
(dofs_per_cell);
//std::vector<Tensor<1,dim> > phi_u (dofs_per_cell);
//std::vector<Tensor<1,dim> > grad_phi_p (dofs_per_cell);
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
scratch.darcy_fe_values.reinit (cell);
typename DoFHandler<dim>::active_cell_iterator
saturation_cell (&triangulation,
cell->level(),
cell->index(),
&saturation_dof_handler);
scratch.saturation_fe_values.reinit (saturation_cell);
//for (; cell!=endc; ++cell, ++saturation_cell)
//{
//darcy_fe_values.reinit (cell);
//saturation_fe_values.reinit (saturation_cell);
data.local_matrix = 0;
scratch.saturation_fe_values.get_function_values
(saturation_solution, scratch.saturation_values);
k_inverse.value_list (scratch.darcy_fe_values.get_quadrature_points(),
scratch.k_inverse_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
const double sat_q = scratch.saturation_values[q];
const double inverse_mobility =
mobility_inverse(sat_q,viscosity);
const double mobility = 1.0 / inverse_mobility;
const Tensor<2,dim> permeability =
invert(scratch.k_inverse_values[q]);
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
scratch.phi_u[k] =
scratch.darcy_fe_values[velocities].value (k,q);
scratch.grad_phi_p[k] =
scratch.darcy_fe_values[pressure].gradient (k,q);
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
data.local_matrix(i,j) += (scratch.k_inverse_values[q] *
inverse_mobility *
scratch.phi_u[i] *
scratch.phi_u[j]
+
permeability * mobility *
scratch.grad_phi_p[i] *
scratch.grad_phi_p[j])
* scratch.darcy_fe_values.JxW(q);
}
}
cell->get_dof_indices (data.local_dof_indices);
//darcy_preconditioner_constraints.distribute_local_to_global
(local_matrix,
// local_dof_indices, darcy_preconditioner_matrix);
//}
}
template <int dim>
void
TwoPhaseFlowProblem<dim>::
copy_local_to_global_darcy_preconditioner (const
Assembly::CopyData::DarcyPreconditioner<dim> &data)
{
darcy_preconditioner_constraints.distribute_local_to_global
(data.local_matrix,
data.local_dof_indices,
darcy_preconditioner_matrix);
}
template <int dim>
void
TwoPhaseFlowProblem<dim>::assemble_darcy_preconditioner ()
{
darcy_preconditioner_matrix = 0;
//const QGauss<dim>
quadrature_formula(parameters.stokes_velocity_degree+1);
const QGauss<dim> quadrature_formula(darcy_degree+2);
typedef
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
CellFilter;
WorkStream::
run (CellFilter (IteratorFilters::LocallyOwnedCell(),
darcy_dof_handler.begin_active()),
CellFilter (IteratorFilters::LocallyOwnedCell(),
darcy_dof_handler.end()),
std_cxx11::bind (&TwoPhaseFlowProblem<dim>::
local_assemble_darcy_preconditioner,
this,
std_cxx11::_1,
std_cxx11::_2,
std_cxx11::_3),
std_cxx11::bind (&TwoPhaseFlowProblem<dim>::
copy_local_to_global_darcy_preconditioner,
this,
std_cxx11::_1),
Assembly::Scratch::
DarcyPreconditioner<dim> (darcy_fe, quadrature_formula,
(update_JxW_values|
update_values|
update_gradients|
update_quadrature_points),
saturation_fe,
update_values),
Assembly::CopyData::
DarcyPreconditioner<dim> (darcy_fe));
darcy_preconditioner_matrix.compress(VectorOperation::add);
}
template <int dim>
void
TwoPhaseFlowProblem<dim>::build_darcy_preconditioner ()
{
computing_timer.enter_section (" Build Darcy preconditioner");
pcout << " Rebuilding Darcy preconditioner..." << std::flush;
assemble_darcy_preconditioner ();
Amg_preconditioner =
std_cxx11::shared_ptr<TrilinosWrappers::PreconditionIC>
(new TrilinosWrappers::PreconditionIC());
Amg_preconditioner->initialize(darcy_preconditioner_matrix.block(0,0));
Mp_preconditioner =
std_cxx11::shared_ptr<TrilinosWrappers::PreconditionIC>
(new TrilinosWrappers::PreconditionIC());
Mp_preconditioner->initialize(darcy_preconditioner_matrix.block(1,1));
pcout << std::endl;
computing_timer.exit_section();
}
template <int dim>
void TwoPhaseFlowProblem<dim>::
local_assemble_darcy_system (const typename
DoFHandler<dim>::active_cell_iterator &cell,
Assembly::Scratch::DarcySystem<dim>
&scratch,
Assembly::CopyData::DarcySystem<dim> &data)
{ //darcy_matrix = 0;
//darcy_rhs = 0;
/*QGauss<dim> quadrature_formula(darcy_degree+2);
QGauss<dim-1> face_quadrature_formula(darcy_degree+2);
FEValues<dim> darcy_fe_values (darcy_fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points |
update_JxW_values);
FEValues<dim> saturation_fe_values (saturation_fe, quadrature_formula,
update_values);
FEFaceValues<dim> darcy_fe_face_values (darcy_fe,
face_quadrature_formula,
update_values |
update_normal_vectors |
update_quadrature_points |
update_JxW_values);*/
//scratch.stokes_fe_values.get_fe().dofs_per_cell;
const unsigned int dofs_per_cell =
scratch.darcy_fe_values.get_fe().dofs_per_cell;
const unsigned int n_q_points =
scratch.darcy_fe_values.n_quadrature_points;
const unsigned int n_face_q_points =
scratch.darcy_fe_face_values.n_quadrature_points;
//FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
//Vector<double> local_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
const PressureRightHandSide<dim> pressure_right_hand_side;
const PressureBoundaryValues<dim> pressure_boundary_values;
/* std::vector<double> pressure_rhs_values (n_q_points);
std::vector<double> boundary_values (n_face_q_points);
std::vector<Tensor<2,dim> > k_inverse_values (n_q_points);
std::vector<double> saturation_values (n_q_points);
std::vector<Tensor<1,dim> > phi_u (dofs_per_cell);
std::vector<double> div_phi_u (dofs_per_cell);
std::vector<double> phi_p (dofs_per_cell);*/
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
scratch.darcy_fe_values.reinit (cell);
typename DoFHandler<dim>::active_cell_iterator
saturation_cell (&triangulation,
cell->level(),
cell->index(),
&saturation_dof_handler);
scratch.saturation_fe_values.reinit (saturation_cell);
//if (rebuild_stokes_matrix)
data.local_matrix = 0;
data.local_rhs = 0;
/*typename DoFHandler<dim>::active_cell_iterator
cell = darcy_dof_handler.begin_active(),
endc = darcy_dof_handler.end();
typename DoFHandler<dim>::active_cell_iterator
saturation_cell = saturation_dof_handler.begin_active();*/
//for (; cell!=endc; ++cell, ++saturation_cell)
//{
//darcy_fe_values.reinit (cell);
//saturation_fe_values.reinit (saturation_cell);
//local_matrix = 0;
//local_rhs = 0;
scratch.saturation_fe_values.get_function_values (saturation_solution,
scratch.saturation_values);
pressure_right_hand_side.value_list
(scratch.darcy_fe_values.get_quadrature_points(),
scratch.pressure_rhs_values);
k_inverse.value_list (scratch.darcy_fe_values.get_quadrature_points(),
scratch.k_inverse_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
scratch.phi_u[k] = scratch.darcy_fe_values[velocities].value
(k,q);
scratch.div_phi_u[k] =
scratch.darcy_fe_values[velocities].divergence (k,q);
scratch.phi_p[k] = scratch.darcy_fe_values[pressure].value
(k,q);
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const double sat_q = scratch.saturation_values[q];
for (unsigned int j=0; j<=i; ++j)
{
data.local_matrix(i,j) += (scratch.phi_u[i] *
scratch.k_inverse_values[q] *
mobility_inverse(sat_q,viscosity)
* scratch.phi_u[j]
- scratch.div_phi_u[i] *
scratch.phi_p[j]
- scratch.phi_p[i] *
scratch.div_phi_u[j])
* scratch.darcy_fe_values.JxW(q);
}
data.local_rhs(i) += (-scratch.phi_p[i] *
scratch.pressure_rhs_values[q])*
scratch.darcy_fe_values.JxW(q);
}
}
for (unsigned int face_no=0;
face_no<GeometryInfo<dim>::faces_per_cell;
++face_no)
if (cell->at_boundary(face_no))
{
scratch.darcy_fe_face_values.reinit (cell, face_no);
pressure_boundary_values
.value_list (scratch.darcy_fe_face_values.get_quadrature_points(),
scratch.boundary_values);
for (unsigned int q=0; q<n_face_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
//const Tensor<1,dim>
scratch.phi_i_u[i] =
scratch.darcy_fe_face_values[velocities].value (i, q);
data.local_rhs(i) += -(scratch.phi_i_u[i] *
scratch.darcy_fe_face_values.normal_vector(q) *
scratch.boundary_values[q] *
scratch.darcy_fe_face_values.JxW(q));
}
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=i+1; j<dofs_per_cell; ++j)
data.local_matrix(i,j) = data.local_matrix(j,i);
cell->get_dof_indices (data.local_dof_indices);
/*darcy_constraints.distribute_local_to_global (local_matrix,
local_rhs,
local_dof_indices,
darcy_matrix,
darcy_rhs);*/
//}
}
template <int dim>
void
TwoPhaseFlowProblem<dim>::
copy_local_to_global_darcy_system (const
Assembly::CopyData::DarcySystem<dim> &data)
{
darcy_constraints.distribute_local_to_global (data.local_matrix,
data.local_rhs,
data.local_dof_indices,
darcy_matrix,
darcy_rhs);
}
template <int dim>
void TwoPhaseFlowProblem<dim>::assemble_darcy_system ()
{
computing_timer.enter_section (" Assemble Darcy system");
darcy_matrix=0;
darcy_rhs=0;
QGauss<dim> quadrature_formula(darcy_degree+2);
QGauss<dim-1> face_quadrature_formula(darcy_degree+2);
//const QGauss<dim>
quadrature_formula(parameters.stokes_velocity_degree+1);
typedef
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
CellFilter;
WorkStream::
run (CellFilter (IteratorFilters::LocallyOwnedCell(),
darcy_dof_handler.begin_active()),
CellFilter (IteratorFilters::LocallyOwnedCell(),
darcy_dof_handler.end()),
std_cxx11::bind (&TwoPhaseFlowProblem<dim>::
local_assemble_darcy_system,
this,
std_cxx11::_1,
std_cxx11::_2,
std_cxx11::_3),
std_cxx11::bind (&TwoPhaseFlowProblem<dim>::
copy_local_to_global_darcy_system,
this,
std_cxx11::_1),
Assembly::Scratch::
DarcySystem<dim> (darcy_fe, quadrature_formula,
face_quadrature_formula,
(update_values|update_normal_vectors|
update_quadrature_points|
update_JxW_values|
update_gradients),
saturation_fe,
update_values),
Assembly::CopyData::
DarcySystem<dim> (darcy_fe));
//if (rebuild_stokes_matrix == true)
darcy_matrix.compress(VectorOperation::add);
darcy_rhs.compress(VectorOperation::add);
//rebuild_stokes_matrix = false;
pcout << std::endl;
computing_timer.exit_section();
std::cout << "assemble ok" << std::endl;
}
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.
