I am working on my graduation thesis using deal.ii, coupling geomechanics
and fluid model based on step-21.
FESystem type is :
fe (FE_RaviartThomas<dim>(degree), 1,
FE_Q<dim>(degree), dim,
FE_DGQ<dim>(degree), 1,
FE_Q<dim>(degree), 1
)
in domain [-1,1]^(dim).
I use normal refinement method so these are no hanging points.
After assemble the system_matrix, I found that I can't use CG method to solve
this problem.
However, When I use UMFPACK's direct solver,
the system shows that :
An error occurred in line <293> of file
</home/chenxi/Downloads/dealii-9.0.1/source/lac/sparse_direct.cc> in function
void dealii::SparseDirectUMFPACK::factorize(const Matrix&) [with Matrix =
dealii::SparseMatrix<double>]
The violated condition was:
status == UMFPACK_OK
Additional information:
UMFPACK routine umfpack_dl_numeric returned error status 1.
I have tried that if I only apply some random number which not equal to zero
to diagnal position, then the solver will work well, however, if I put number
in other places,
this will show an error.
Thank you in advance.
(PS: I am a newer to Linux and starting trying Deal.II. I have checked the code
in sparse_direct.cc sparse_direct.cc and umfpack.h, but I still can't find out
the mistake)
--
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.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/base/symmetric_tensor.h>
#include <fstream>
#include <sstream>
#include <iostream>
#include <fstream>
#include <deal.II/base/tensor_function.h>
namespace Step21
{
using namespace dealii;
template <int dim>
class TwoPhaseFlowProblem
{
public:
TwoPhaseFlowProblem (const unsigned int degree);
void run ();
private:
void make_grid_and_dofs ();
void assemble_system ();
void assemble_rhs_S ();
double get_maximal_velocity () const;
void solve ();
void project_back_saturation ();
void output_results () const;
const unsigned int degree;
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern; //update-code-8
//DynamicSparsityPattern dsp;
SparseMatrix<double> system_matrix;
ConstraintMatrix constraints;
const unsigned int n_refinement_steps;
double time_step;
unsigned int timestep_number;
double viscosity;
Vector<double> solution;
Vector<double> old_solution;
Vector<double> system_rhs;
};
template <int dim>
class PressureRightHandSide : public Function<dim>
{
public:
PressureRightHandSide () : Function<dim>(1) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double
PressureRightHandSide<dim>::value (const Point<dim> &/*p*/,
const unsigned int /*component*/) const
{
return 0;
}
template <int dim>
class PressureBoundaryValues : public Function<dim>
{
public:
PressureBoundaryValues () : Function<dim>(1) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double
PressureBoundaryValues<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
return 1-p[0];
}
template <int dim>
class SaturationBoundaryValues : public Function<dim>
{
public:
SaturationBoundaryValues () : Function<dim>(1) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double
SaturationBoundaryValues<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
if (p[0] == 0)
return 1;
else
return 0;
}
template <int dim> //update-code
class DisplacementBoundaryValues : public Function<dim> //update-code
{ //update-code
public:
DisplacementBoundaryValues () : Function<dim>(dim) {} //update-code
//update-code
virtual double value (const Point<dim> &p, //update-code
const unsigned int component) const; //update-code
virtual void vector_value (const Point<dim> &p, //update-code
Vector<double> &value) const; //update-code
}; //update-code
template<int dim> //update-code
double //update-code
DisplacementBoundaryValues<dim>::value(const Point<dim> &p, //update-code
const unsigned int component) const //update-code
{ //update-code
return 0; //update-code
} //update-code
template<int dim> //update-code
void //update-code
DisplacementBoundaryValues<dim>::vector_value(const Point<dim> &p, //update-code
Vector<double> &value) const //update-code
{ //update-code
ZeroFunction<dim>(dim).vector_value(p,value); //update-code //update-code
} //update-code
template <int dim>
class InitialValues : public Function<dim>
{
public:
InitialValues () : Function<dim>(2*dim+2) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
};
template <int dim>
double
InitialValues<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
return ZeroFunction<dim>(2*dim+2).value (p, component);
}
template <int dim>
void
InitialValues<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
ZeroFunction<dim>(2*dim+2).vector_value (p, values);
}
namespace SingleCurvingCrack
{
template <int dim>
class KInverse : public TensorFunction<2,dim>
{
public:
KInverse ()
:
TensorFunction<2,dim> ()
{}
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<Tensor<2,dim> > &values) const;
};
template <int dim>
void
KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
std::vector<Tensor<2,dim> > &values) const
{
Assert (points.size() == values.size(),
ExcDimensionMismatch (points.size(), values.size()));
for (unsigned int p=0; p<points.size(); ++p)
{
values[p].clear ();
const double distance_to_flowline
= std::fabs(points[p][1]-0.5-0.1*std::sin(10*points[p][0]));
const double permeability = std::max(std::exp(-(distance_to_flowline*
distance_to_flowline)
/ (0.1 * 0.1)),
0.01);
for (unsigned int d=0; d<dim; ++d)
values[p][d][d] = 1./permeability;
}
}
}
namespace RandomMedium
{
template <int dim>
class KInverse : public TensorFunction<2,dim>
{
public:
KInverse ()
:
TensorFunction<2,dim> ()
{}
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<Tensor<2,dim> > &values) const;
private:
static std::vector<Point<dim> > centers;
static std::vector<Point<dim> > get_centers ();
};
template <int dim>
std::vector<Point<dim> >
KInverse<dim>::centers = KInverse<dim>::get_centers();
template <int dim>
std::vector<Point<dim> >
KInverse<dim>::get_centers ()
{
const unsigned int N = (dim == 2 ?
40 :
(dim == 3 ?
100 :
throw ExcNotImplemented()));
std::vector<Point<dim> > centers_list (N);
for (unsigned int i=0; i<N; ++i)
for (unsigned int d=0; d<dim; ++d)
centers_list[i][d] = static_cast<double>(rand())/RAND_MAX;
return centers_list;
}
template <int dim>
void
KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
std::vector<Tensor<2,dim> > &values) const
{
Assert (points.size() == values.size(),
ExcDimensionMismatch (points.size(), values.size()));
for (unsigned int p=0; p<points.size(); ++p)
{
values[p].clear ();
double permeability = 0;
for (unsigned int i=0; i<centers.size(); ++i)
permeability += std::exp(-(points[p]-centers[i])*(points[p]-centers[i])
/ (0.05 * 0.05));
const double normalized_permeability
= std::min (std::max(permeability, 0.01), 4.);
for (unsigned int d=0; d<dim; ++d)
values[p][d][d] = 1./normalized_permeability;
}
}
}
double mobility_inverse (const double S,
const double viscosity)
{
return 1.0 / (1.0/viscosity * S * S + (1-S) * (1-S));
}
double fractional_flow (const double S,
const double viscosity)
{
return S*S / (S * S + viscosity * (1-S) * (1-S));
}
template <class Matrix>
class InverseMatrix : public Subscriptor
{
public:
InverseMatrix (const Matrix &m);
void vmult (Vector<double> &dst,
const Vector<double> &src) const;
private:
const SmartPointer<const Matrix> matrix;
};
template <class Matrix>
InverseMatrix<Matrix>::InverseMatrix (const Matrix &m)
:
matrix (&m)
{}
template <class Matrix>
void InverseMatrix<Matrix>::vmult (Vector<double> &dst,
const Vector<double> &src) const
{
SolverControl solver_control (std::max<unsigned int>(src.size(), 200),
1e-8*src.l2_norm());
SolverCG<> cg (solver_control);
dst = 0;
cg.solve (*matrix, dst, src, PreconditionIdentity());
}
template <int dim>
TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem (const unsigned int degree)
:
degree (degree),
fe (FE_RaviartThomas<dim>(degree), 1,
FE_Q<dim>(degree), dim,
FE_DGQ<dim>(degree), 1,
FE_Q<dim>(degree), 1
), // update-code-1
dof_handler (triangulation),
n_refinement_steps (4),
time_step (0),
viscosity (0.2)
{}
template <int dim>
void TwoPhaseFlowProblem<dim>::make_grid_and_dofs ()
{
GridGenerator::hyper_cube (triangulation, 0, 1);
triangulation.refine_global (n_refinement_steps);
dof_handler.distribute_dofs (fe);
constraints.clear();
//DoFRenumbering::component_wise (dof_handler);
//std::vector<types::global_dof_index> dofs_per_component (2*dim+2);
//DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);
//const unsigned int n_u = dofs_per_component[0],
// n_displacement = dofs_per_component[dim], //update-code-5
// n_p = dofs_per_component[2*dim],
// n_s = dofs_per_component[2*dim+1];
DoFTools::make_hanging_node_constraints(dof_handler,constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs(),dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ true);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
old_solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
template <int dim>
void TwoPhaseFlowProblem<dim>::assemble_system ()
{
system_matrix=0;
system_rhs=0;
QGauss<dim> quadrature_formula(degree+2);
QGauss<dim-1> face_quadrature_formula(degree+2);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
update_values | update_normal_vectors |
update_quadrature_points | update_JxW_values);
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 = face_quadrature_formula.size();
const double lambda = 1.0; //update-code-7
const double biot = 0.8;
const double phi = 0.2;
const double K_s = 0.2;
const double c_w = 0.9;
const double c_o = 1.1;
Tensor<1,dim> t;
for(int ii=0;ii<dim;++ii)
t[ii]=100000;
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;
const RandomMedium::KInverse<dim> k_inverse;
const DisplacementBoundaryValues<dim> displacement_boundary_values; //update-code-6
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<Vector<double> > old_solution_values(n_q_points, Vector<double>(2*dim+2));
std::vector<SymmetricTensor<2,dim> > old_solution_values_symmetricgradient_displacement(n_q_points);
std::vector<double> old_solution_values_divergence_displacement(n_q_points);
std::vector<double> old_solution_values_divergence_velocities(n_q_points);
std::vector<Tensor<1,dim> > old_solution_values_gradient_saturation(n_q_points);
const FEValuesExtractors::Vector velocities (0); //update-code-5
const FEValuesExtractors::Vector displacement(dim);
const FEValuesExtractors::Scalar pressure (2*dim);
const FEValuesExtractors::Scalar saturation (2*dim+1);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
local_matrix = 0;
local_rhs = 0;
fe_values.get_function_values (old_solution, old_solution_values);
fe_values[displacement].get_function_divergences(old_solution,old_solution_values_divergence_displacement);
fe_values[displacement].get_function_symmetric_gradients(old_solution,old_solution_values_symmetricgradient_displacement);
fe_values[velocities].get_function_divergences(old_solution,old_solution_values_divergence_velocities);
fe_values[pressure].get_function_gradients(old_solution,old_solution_values_gradient_saturation);
fe_values[saturation].get_function_gradients(old_solution,old_solution_values_gradient_saturation);
pressure_right_hand_side.value_list (fe_values.get_quadrature_points(),
pressure_rhs_values);
k_inverse.value_list (fe_values.get_quadrature_points(),
k_inverse_values);
//sigma zero
Tensor<2,dim> sigma_zero;
for(int ii=0;ii<dim;++ii)
{
for(int kk = 0;kk<dim;++kk)
{
sigma_zero[ii][kk] = 10000;
}
}
//p-zero
const double p_zero = 100000;
for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const double old_p = old_solution_values[q](2*dim);
//const Tensor<1,dim> old_displacement;
//for(int ii=0;ii<dim;++ii)
// old_displacement[ii] = old_solution_values[q][dim+ii];
const double old_s = old_solution_values[q](2*dim+1);
const double N = (biot-phi)/K_s + old_s*phi*c_w + (1-old_s)*phi*c_o;
const double M = old_s*((biot-phi)/K_s+phi*c_w);
const double Fw = fractional_flow(old_s,viscosity);
Tensor<1,dim> t;
for(int ii = 0;ii<dim;++ii)
t[ii] = 10000;
Tensor<1,dim> old_u;
for(int ii=0;ii<dim;++ii)
old_u[ii]=old_solution_values[q][ii];
const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q);
const Tensor<1,dim> phi_i_displacement = fe_values[displacement].value (i, q);
const Tensor<2,dim> sigma_phi_i_displacement = fe_values[displacement].gradient(i,q);
const double div_phi_i_u = fe_values[velocities].divergence (i, q);
const double div_phi_i_displacement = fe_values[displacement].divergence(i,q);
//const double div_phi_i_p = fe_values[pressure].divergence(i,q);
const double phi_i_p = fe_values[pressure].value (i, q);
const double phi_i_s = fe_values[saturation].value (i, q);
const Tensor<1,dim> grad_phi_i_s = fe_values[saturation].gradient(i,q);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const Tensor<1,dim> phi_j_u = fe_values[velocities].value (j, q);
const Tensor<1,dim> phi_j_displacement = fe_values[displacement].value (j, q);
const Tensor<2,dim> sigma_phi_j_displacement = fe_values[displacement].gradient(j,q);
const double div_phi_j_displacement = fe_values[displacement].divergence(j,q);
const double div_phi_j_u = fe_values[velocities].divergence (j, q);
const double phi_j_p = fe_values[pressure].value (j, q);
const double phi_j_s = fe_values[saturation].value (j, q);
const Tensor<1,dim> grad_phi_j_p = fe_values[pressure].gradient(j,q);
const Tensor<1,dim> grad_phi_j_s = fe_values[saturation].gradient(j,q);
//if(i==j) if apply this condition then the solver will work well!
local_matrix(i,j) += (1)
* fe_values.JxW(q);
}
local_rhs(i) += (1)*
fe_values.JxW(q);
}
//std::cout<<std::endl<<std::endl;
for (unsigned int face_no=0;
face_no<GeometryInfo<dim>::faces_per_cell;
++face_no)
if (cell->at_boundary(face_no))
{
fe_face_values.reinit (cell, face_no);
pressure_boundary_values
.value_list (fe_face_values.get_quadrature_points(),
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>
phi_i_displacement = fe_face_values[displacement].value (i, q);
local_rhs(i) += -(1)
* fe_face_values.JxW(q);
}
//local_rhs(i)+=6;
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
local_matrix(i,j));
for (unsigned int i=0; i<dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += local_rhs(i);
}
}
template <int dim>
void TwoPhaseFlowProblem<dim>::solve ()
{
SparseDirectUMFPACK direct_solver;
direct_solver.initialize(system_matrix);
direct_solver.vmult(solution,system_rhs);
//constriants.distribute(solution);
// Finally, we have to take care of the saturation equation. The first
// business we have here is to determine the time step using the formula
// in the introduction. Knowing the shape of our domain and that we
// created the mesh by regular subdivision of cells, we can compute the
// diameter of each of our cells quite easily (in fact we use the linear
// extensions in coordinate directions of the cells, not the
// diameter). Note that we will learn a more general way to do this in
// step-24, where we use the GridTools::minimal_cell_diameter function.
//
// The maximal velocity we compute using a helper function to compute the
// maximal velocity defined below, and with all this we can evaluate our
// new time step length:
time_step = 0.01;
// The next step is to assemble the right hand side, and then to pass
// everything on for solution. At the end, we project back saturations
// onto the physically reasonable range:
old_solution = solution;
}
template <int dim>
void TwoPhaseFlowProblem<dim>::output_results () const
{
if (timestep_number % 5 != 0)
return;
std::vector<std::string> solution_names;
switch (dim)
{
case 2:
solution_names.push_back ("u");
solution_names.push_back ("v");
solution_names.push_back ("displacement_u");
solution_names.push_back ("displacement_v");
solution_names.push_back ("p");
solution_names.push_back ("S");
break;
case 3:
solution_names.push_back ("u");
solution_names.push_back ("v");
solution_names.push_back ("w");
solution_names.push_back ("displacement_u");
solution_names.push_back ("displacement_v");
solution_names.push_back ("displacement_w");
solution_names.push_back ("p");
solution_names.push_back ("S");
break;
default:
Assert (false, ExcNotImplemented());
}
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, solution_names);
data_out.build_patches (degree+1);
std::ostringstream filename;
filename << "solution-" << timestep_number << ".vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
template <int dim>
void TwoPhaseFlowProblem<dim>::run ()
{
make_grid_and_dofs();
{
ConstraintMatrix constraints;
constraints.close();
VectorTools::project (dof_handler,
constraints,
QGauss<dim>(degree+2),
InitialValues<dim>(),
old_solution);
}
timestep_number = 1;
double time = 0;
do
{
std::cout << "Timestep " << timestep_number
<< std::endl;
assemble_system ();
solve ();
output_results ();
time += time_step;
++timestep_number;
std::cout << " Now at t=" << time
<< ", dt=" << time_step << '.'
<< std::endl
<< std::endl;
}
while (time <= 1.);
}
}
int main ()
{
try
{
using namespace dealii;
using namespace Step21;
deallog.depth_console (0);
TwoPhaseFlowProblem<2> two_phase_flow_problem(1);
two_phase_flow_problem.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
