[Trilinos-Users] Two solvers (direct and iterative) to solve the matrices created through Newton method
Isabel Gil
misabel at litec.csic.es
Fri Sep 12 05:26:54 MDT 2008
Hi Daniele and everybody,
first of all, thanks for your reply. I write my doubts in the middle of
your previous email just in case you could help me.
Really I am using a Finite Element library and together with this one I
use Trilinos because I want to use automatic differentiation that
provides the package Sacado.
(Between //////////////////////// there are the lines of code that have
been programed)
Daniele Avitabile wrote:
> Isa,
>
> if I were you I would make sure that you are solving the right linear
> system. The fact that your direct solver doesn't converge could mean
> that something is wrong in the outer loop (i.e. the Newton loop). Pay
> special attention to the way you pass x and b... what are x and b in
> your code? Just print them out, together with the matrix. You can have
> them printed on the screen via cout.
I write some parts of the code, because I don't know where the problem
could be. But, firstly, I will explain the problem I want to solve:
I want to solve NavierStokes equations and the momentum equations have
to be solved through Newton method. For that I am using automatic
differentiation package, Sacado.
I will use Newton method ONLY for momentum equations. For pressure
equation I will use another method. I mean, the residual matrix to use
Newton method only contains velocity nodes, not pressure ones (well,
that is what I want to get.)
I am using a finite element library: I use Q2 finite element for
velocity and Q1 finite element for pressure, then I create a finite
element system that keep all together. DIMENSION=2 because we work with
(u,v) as velocity vector:
////////////////////////////////////////////////////////////////////////////////
fe_velocity(2)
fe_pressure(1)
fe_total(fe_velocity,DIMENSION, fe_pressure,1)
std::auto_ptr<Epetra_Map> Map;
std::auto_ptr<Epetra_CrsMatrix> Matrix;
Vector<double> right_hand_side;
Vector<double> old_solution;
Vector<double> current_solution;
Vector<double> predictor;
////////////////////////////////////////////////////////////////////////////////
It is supposed that "Matrix" and "right_hand_side" are only related to
velocity nodes (nodes of component u + nodes of component v), NOT to
total nodes (nodes of component u + nodes of component v + nodes of
pressure p)
I want to create a "Matrix" that only has the derivatives of the
residual in relation to velocity components. But I don't know if the
code below is
the correct one:
First I create the "Epetra Matrix":
////////////////////////////////////////////////////////////////////////////////
void setup_system()
{
//size[0]=the number of velocity nodes.
Map.reset (new Epetra_Map(size[0], 0, communicator));
std::vector<int> row_lengths (size[0]);
//sp_general.block(0,0) is a sparsity pattern that only has into count
velocity nodes.
for (unsigned int i=0; i<size[0]; ++i)
row_lengths[i] = sp_general.block(0,0).row_length (i);
Matrix.reset (new Epetra_CrsMatrix(Copy, *Map, &row_lengths[0], true));
const unsigned int max_nonzero_entries
= *std::max_element (row_lengths.begin(), row_lengths.end());
std::vector<double> values(max_nonzero_entries, 0);
std::vector<int> row_indices(max_nonzero_entries);
for (unsigned int row=0; row<size[0]; ++row)
{
row_indices.resize (row_lengths[row], 0);
values.resize (row_lengths[row], 0.);
for (int i=0; i<row_lengths[row]; ++i)
row_indices[i] = sp_general.block(0,0).column_number (row, i);
Matrix->InsertGlobalValues(row, row_lengths[row],
&values[0], &row_indices[0]);
}
Matrix->FillComplete();
}
////////////////////////////////////////////////////////////////////////////////
Secondly, I create this matrix:
////////////////////////////////////////////////////////////////////////////////
void assemble_system ()
{
// nodes of a cell= nodes of velocity + nodes of pressure.
const unsigned int dofs_per_cell =
dof_handler_total.get_fe().dofs_per_cell;
// nodes of a cell= nodes of velocity, for pressure we are using Q1, so
nodes of pressure=4
const unsigned int velocity_dofs_per_cell=dofs_per_cell-4;
std::vector<unsigned int> dof_indices (dofs_per_cell);
std::vector<unsigned int> dof_indices_neighbor (dofs_per_cell);
const UpdateFlags update_flags = update_values
| update_gradients
| update_q_points
| update_JxW_values;
FEValues<DIMENSION> fe_v (mapping, fe_total,
quadrature,
update_flags);
DoFHandler<DIMENSION>::active_cell_iterator
cell = dof_handler_total.begin_active(),
endc = dof_handler_total.end();
for (; cell!=endc; ++cell)
{
fe_v.reinit (cell);
cell->get_dof_indices (dof_indices);
assemble_cell_term(fe_v, dof_indices);
}
Matrix->FillComplete();
}
////////////////////////////////////////////////////////////////////////////////
void assemble_cell_term (const FEValues<DIMENSION> &fe_v,
const std::vector<unsigned int> &dof_indices)
{
//All nodes: the two components of velocity (u,v) and the pressure ones.
const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
const unsigned int n_q_points = fe_v.n_quadrature_points;
//Only velocity nodes.
const unsigned int velocity_dofs_per_cell=dofs_per_cell;
Table<2,Sacado::Fad::DFad<double> >
W (n_q_points, DIMENSION+1);
Table<2,double>
W_old (n_q_points, DIMENSION+1);
Table<2,Sacado::Fad::DFad<double> >
W_theta (n_q_points, DIMENSION+1);
Table<3,Sacado::Fad::DFad<double> >
grad_W (n_q_points, DIMENSION+1, DIMENSION);
std::vector<Sacado::Fad::DFad<double> >
independent_local_dof_values(dofs_per_cell);
for (unsigned int i=0; i<dofs_per_cell; ++i)
independent_local_dof_values[i] = current_solution(dof_indices[i]);
for (unsigned int i=0; i<dofs_per_cell; ++i)
independent_local_dof_values[i].diff(i, dofs_per_cell);
for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int c=0; c<DIMENSION+1; ++c)
{
W[q][c] = 0;
W_old[q][c] = 0;
W_theta[q][c] = 0;
for (unsigned int d=0; d<DIMENSION; ++d)
grad_W[q][c][d] = 0;
}
for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int c =
fe_v.get_fe().system_to_component_index(i).first;
W[q][c] += independent_local_dof_values[i] *
fe_v.shape_value_component(i, q, c);
W_old[q][c] += old_solution(dof_indices[i]) *
fe_v.shape_value_component(i, q, c);
W_theta[q][c] += (theta *
independent_local_dof_values[i]
+
(1-theta) *
old_solution(dof_indices[i])) *
fe_v.shape_value_component(i, q, c);
for (unsigned int d = 0; d < DIMENSION; d++)
grad_W[q][c][d] += independent_local_dof_values[i] *
fe_v.shape_grad_component(i, q, c)[d];
}
std::vector<double> residual_derivatives (velocity_dofs_per_cell);
for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
{
Sacado::Fad::DFad<double> F_i = 0;
const unsigned int
component_i = fe_v.get_fe().system_to_component_index(i).first;
if (component_i!=DIMENSION)
{
for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
{
//Temporal part:
F_i += ((1.0 /time_step_value) *
(W[point][component_i] - W_old[point][component_i]) *
fe_v.shape_value_component(i, point, component_i)*
fe_v.JxW(point));
//Gradient of pressure:
F_i-=(fe_v.shape_grad_component(i, point,
component_i)[component_i] *
W_old[point][DIMENSION]*
fe_v.JxW(point)
);
//Diffusion part:
for (unsigned int d=0; d<DIMENSION; d++)
F_i+=(fe_v.JxW(point)*viscosity*fe_v.shape_grad_component(i, point,
component_i)[d] *
grad_W[point][component_i][d]
);
//Convection part:
for (unsigned int k=0;k<DIMENSION;k++)
F_i+=
(fe_v.JxW(point)*W[point][k]*grad_W[point][component_i][k]*
fe_v.shape_value_component(i, point, component_i));
}
for (unsigned int k=0; k<velocity_dofs_per_cell; ++k)
residual_derivatives[k] = F_i.fastAccessDx(k);
Matrix->SumIntoGlobalValues(dof_indices[i],
velocity_dofs_per_cell,
&residual_derivatives[0],
reinterpret_cast<int*>(
const_cast<unsigned int*>(
&dof_indices[0])));
right_hand_side(dof_indices[i]) -= F_i.val();
}
}
}
////////////////////////////////////////////////////////////////////////////////
Here the outer loop comes:
////////////////////////////////////////////////////////////////////////////////
while (iteration <= num_iterations)
{
unsigned int nonlin_iter = 0;
current_solution = predictor;
while (true)
{
Matrix->PutScalar(0);
Matrix->FillComplete();
right_hand_side = 0;
assemble_system ();
const double res_norm = right_hand_side.l2_norm();
if (std::fabs(res_norm) < 1e-10)
{
std::printf(" %-16.3e (converged)\n\n", res_norm);
break;
}
else
{
newton_update = 0;
std::pair<unsigned int, double> convergence
= solve (newton_update);
for (unsigned int i=0;i<size[0];i++)
current_solution(i) += newton_update(i);
std::printf(" %-16.3e %04d %-5.2e\n",res_norm,
convergence.first, convergence.second);
}
++nonlin_iter;
AssertThrow (nonlin_iter <= 10,
ExcMessage ("No convergence in nonlinear solver"));
}
predictor = current_solution;
predictor.sadd (2.0, -1.0, old_solution);
old_solution = current_solution;
///////////////////////////////////////////////////////////////////////////////
I want to create a "Matrix" that only has the derivatives of the
residual in relation to velocity components. But, obviously, I am doing
something wrong when I create the matrix.
the correct
>
> The error in the iterative solver is probably because your comm is
> *not* an Epetra_MPIComm, as you define it to be an Epetra_Comm. Let me
> know if you are not familiar with MPI.
I don't understand very well what you want to say.
I have done the configure putting --enable-mpi --with-mpi-compilers.
Later I use :
#include <Epetra_SerialComm.h>
Epetra_SerialComm communicator;
std::auto_ptr<Epetra_Map> Map;
std::auto_ptr<Epetra_CrsMatrix> Matrix;
So, I don't know what I should do to use the iterative solver.
Thank-you very much in advance.
Any help will be welcome.
Best regards.
Isa
>
> Best.
>
> Daniele
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://software.sandia.gov/mailman/private/trilinos-users/attachments/20080912/4ec8cb4c/attachment-0001.html
More information about the Trilinos-Users
mailing list