Hi everybody,
I am following the tutorial step-33 in order to code my program.
In my case, I want to solve NavierStokes equations and the momentum
equations have to be solved through Newton iterations. For that I am using
automatic differentiation package, Sacado.
I will use Trilinos algebra only for momentum equations, for pressure
equation I will use the methods that dealII provides.
So, I have:
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 dofs, so, following step-33:
////////////////////////////////////////////////////////////////////////////////
void setup_system ()
{
std::vector< unsigned int> size(2);
std::vector< unsigned int > dofs_per_component (DIMENSION+1);
DoFTools::count_dofs_per_component (dof_handler_total,
dofs_per_component);
size[0]=dofs_per_component[0]+dofs_per_component[1];
size[1]=dofs_per_component[2];
//right_hand_side only works with velocity dofs:
right_hand_side.reinit (size[0]);
//The rest of vectors will save the three components (two velocities and
the pressure)
old_solution.reinit (dof_handler_total.n_dofs());
current_solution.reinit (dof_handler_total.n_dofs());
predictor.reinit (dof_handler_total.n_dofs());
Map.reset (new Epetra_Map(size[0], 0, communicator));
//sp_general.block(0,0) is the sparsity pattern that only works with
velocity dofs.
std::vector<int> row_lengths (size[0]);
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();
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void assemble_system ()
{
const unsigned int dofs_per_cell =
dof_handler_total.get_fe().dofs_per_cell;
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)
{
const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
const unsigned int n_q_points = fe_v.n_quadrature_points;
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 (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;
//I only take the components different from the pressure ones.
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[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<(dofs_per_cell); ++k)
residual_derivatives[k] = F_i.fastAccessDx(k);
Matrix->SumIntoGlobalValues(dof_indices[i],
(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();
}
}
}
////////////////////////////////////////////////////////////////////////////////
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 above code is
the correct
one since when I run the code (obviously there is also a "solve method"
using "Aztec00 solver") I obtain the following problem:
Number of active cells: 1024
Number of degrees of freedom: 9539
NonLin Res Lin Iter Lin Res
_____________________________________
AMESOS ERROR -22, ../../../../packages/amesos/src/Amesos_Klu.cpp, line 545
AMESOS ERROR -22, ../../../../packages/amesos/src/Amesos_Klu.cpp, line 656
AMESOS ERROR -22, ../../../../packages/amesos/src/Amesos_Klu.cpp, line 545
AMESOS ERROR -22, ../../../../packages/amesos/src/Amesos_Klu.cpp, line 656
AMESOS ERROR -22, ../../../../packages/amesos/src/Amesos_Klu.cpp, line 696
2.175e-01 0000 0.00e+00
AMESOS ERROR -22, ../../../../packages/amesos/src/Amesos_Klu.cpp, line 545
AMESOS ERROR -22, ../../../../packages/amesos/src/Amesos_Klu.cpp, line 656
AMESOS ERROR -22, ../../../../packages/amesos/src/Amesos_Klu.cpp, line 545
AMESOS ERROR -22, ../../../../packages/amesos/src/Amesos_Klu.cpp, line 656
AMESOS ERROR -22, ../../../../packages/amesos/src/Amesos_Klu.cpp, line 696
2.175e-01 0000 0.00e+00
.
.
.
I think I am creating a matrix and a right hand side that are not the
actual ones I am looking for since I want the Jacobian matrix of the
residual with respect to velocity vector.
Thanks in advance.
Best
Isa
_______________________________________________
