Could get it to work, now the code works equally fine for one or several
MPI threads. Still, the next open question is how to apply non-zero
Dirichlet boundary conditions (as in the file, for cos(pi*x)*cos(pi*y).
When applying them as I do in the file, I get the warning that the solver
is not converged. When using Neumann-Conditions, the solver converges and I
get the correct output. What would be the solution for that?
Thanks!
Am Montag, 29. April 2019 10:20:13 UTC+2 schrieb Maxi Miller:
>
> Now with file...
>
> Am Montag, 29. April 2019 10:19:44 UTC+2 schrieb Maxi Miller:
>>
>> The functions computePreconditioner and computeJacobian are necessary,
>> else I will get compilation problems. Thus I think the current version is
>> the best MWE I can get. Running with one thread gives a "Test passed",
>> using more threads will result in "Convergence failed".
>>
>> Am Samstag, 27. April 2019 05:43:47 UTC+2 schrieb Wolfgang Bangerth:
>>>
>>> On 4/26/19 2:28 PM, 'Maxi Miller' via deal.II User Group wrote:
>>> > I think the current version is the smallest version I can get. It will
>>> work
>>> > for a single thread, but not for two or more threads.
>>>
>>> I tried to run this on my machine, but I don't have NOX as part of my
>>> Trilinos
>>> installation. That might have to wait for next week.
>>>
>>> But I think there is still room for making the program shorter. The
>>> point is
>>> that the parts of the program that run before the error happens do not
>>> actually have to make sense. For example, could you replace your own
>>> boundary
>>> values class by ZeroFunction? Instead of assembling something real,
>>> could you
>>> just use the identity matrix on each cell?
>>>
>>> Other things you don't need: Timers, convergence tables, etc. Which of
>>> the
>>> parameters you set in solve_for_NOX() are necessary? (Necessary to
>>> reproduce
>>> the bug; I know they're necessary for your application, but that's
>>> irrelevant
>>> here.) Could you replace building the residual by just putting random
>>> stuff in
>>> the vector?
>>>
>>> Similarly, the computeJacobian and following functions really just check
>>> conditions, but when you run the program to illustrate the bug you're
>>> chasing,
>>> do these checks trigger? If not, remove them, then remove the calls to
>>> these
>>> functions (because they don't do anything any more), then remove the
>>> whole
>>> function.
>>>
>>> I think there's still room to make this program smaller by at least a
>>> factor
>>> of 2 or 3 or 4 :-)
>>>
>>> Best
>>> W.
>>>
>>> --
>>> ------------------------------------------------------------------------
>>> Wolfgang Bangerth email: [email protected]
>>> www:
>>> http://www.math.colostate.edu/~bangerth/
>>>
>>>
--
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.
//Nox include files
#include <NOX.H>
#include <NOX_LAPACK_Group.H>
#include <Teuchos_StandardCatchMacros.hpp>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/work_stream.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/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/generic_linear_algebra.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/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_q.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/meshworker/mesh_loop.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/meshworker/scratch_data.h>
#include <deal.II/meshworker/copy_data.h>
// NOX Objects
#include "NOX.H"
#include "NOX_Epetra.H"
// Trilinos Objects
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Vector.h"
#include "Epetra_Operator.h"
#include "Epetra_RowMatrix.h"
#include "NOX_Epetra_Interface_Required.H" // base class
#include "NOX_Epetra_Interface_Jacobian.H" // base class
#include "NOX_Epetra_Interface_Preconditioner.H" // base class
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "Epetra_LinearProblem.h"
#include "AztecOO.h"
#include "Teuchos_StandardCatchMacros.hpp"
// User's application specific files
//#include "problem_interface.h" // Interface file to NOX
//#include "finite_element_problem.h"
// Required for reading and writing parameter lists from xml format
// Configure Trilinos with --enable-teuchos-extended
#ifdef HAVE_TEUCHOS_EXTENDED
#include "Teuchos_XMLParameterListHelpers.hpp"
#include <sstream>
#endif
// This is C++ ...
#include <fstream>
#include <iostream>
using namespace dealii;
template <int dim>
class BoundaryValues : public dealii::Function<dim>
{
public:
BoundaryValues(const size_t n_components)
: dealii::Function<dim>(1), n_components(n_components)
{}
virtual double value (const dealii::Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value(const dealii::Point<dim> &p, dealii::Vector<double> &value) const;
private:
const size_t n_components;
};
template <int dim>
double BoundaryValues<dim>::value (const dealii::Point<dim> &p,
const unsigned int) const
{
const double x = p[0];
const double y = p[1];
return cos(M_PI * x) * cos(M_PI * y);//sin(M_PI * x) * sin(M_PI * y);
}
template <int dim>
void BoundaryValues<dim>::vector_value(const dealii::Point<dim> &p, dealii::Vector<double> &value) const
{
for(size_t i = 0; i < value.size(); ++i)
value[i] = BoundaryValues<dim>::value(p, i);
}
template <int dim>
class ProblemInterface : public NOX::Epetra::Interface::Required,
public NOX::Epetra::Interface::Preconditioner,
public NOX::Epetra::Interface::Jacobian
{
public:
ProblemInterface(std::function<void(const TrilinosWrappers::MPI::Vector&, TrilinosWrappers::MPI::Vector&)> residual_function,
const MPI_Comm &mpi_communicator,
const IndexSet &locally_owned_dofs,
const IndexSet &locally_relevant_dofs) :
residual_function(residual_function),
mpi_communicator(mpi_communicator),
locally_owned_dofs(locally_owned_dofs),
locally_relevant_dofs(locally_relevant_dofs)
{
};
~ProblemInterface(){
};
bool computeF(const Epetra_Vector &x, Epetra_Vector &FVec, NOX::Epetra::Interface::Required::FillType)
{
f_vec.reinit(locally_owned_dofs, mpi_communicator);
residual_vec.reinit(locally_owned_dofs, mpi_communicator);
Epetra_Vector f_epetra_vec = Epetra_Vector(View, f_vec.trilinos_vector(), 0);
residual_vec = 0.;
f_epetra_vec = x;
residual_function(f_vec, residual_vec);
for(auto index : residual_vec.locally_owned_elements())
FVec[locally_owned_dofs.index_within_set(index)] = residual_vec[index];
FVec = Epetra_Vector(Copy, residual_vec.trilinos_vector(), 0);
return true;
}
bool computeJacobian(const Epetra_Vector &x, Epetra_Operator &Jac)
{
Epetra_CrsMatrix* Jacobian = dynamic_cast<Epetra_CrsMatrix*>(&Jac);
if (Jacobian == NULL) {
std::cout << "ERROR: Problem_Interface::computeJacobian() - The supplied"
<< "Epetra_Operator is NOT an Epetra_RowMatrix!" << std::endl;
throw;
}
(void) x;
return false;
}
bool computePreconditioner(const Epetra_Vector& x, Epetra_Operator& M,
Teuchos::ParameterList* precParams = 0)
{
(void) x;
(void) M;
(void) precParams;
std::cout << "ERROR: Problem_Interface::preconditionVector() - Use Explicit Jacobian only for this test problem!" << std::endl;
throw 1;
return false;
}
private:
std::function<void(const TrilinosWrappers::MPI::Vector&, TrilinosWrappers::MPI::Vector&)> residual_function;
Epetra_Vector *rhs;
MPI_Comm mpi_communicator;
IndexSet locally_owned_dofs, locally_relevant_dofs;
TrilinosWrappers::MPI::Vector f_vec, residual_vec;
};
template <int dim>
class Step5
{
public:
Step5();
~Step5();
void run();
private:
void setup_system();
void assemble_system();
void calculate_residual(const TrilinosWrappers::MPI::Vector &u, TrilinosWrappers::MPI::Vector &residual);
void solve_with_NOX();
void assemble_rhs(double &return_value,
const double factor,
const double &val_u,
const Tensor<1, dim> &grad_u,
const Point<dim> p,
const double &fe_values_value,
const Tensor<1, dim> &fe_gradient_value,
const double JxW_value);
MPI_Comm mpi_communicator;
parallel::distributed::Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FE_Q<dim> fe;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
SparsityPattern sparsity_pattern;
TrilinosWrappers::SparseMatrix system_matrix;
TrilinosWrappers::MPI::Vector locally_relevant_solution;
TrilinosWrappers::MPI::Vector system_rhs;
TrilinosWrappers::MPI::Vector evaluation_point;
AffineConstraints<double> hanging_node_constraints;
ConditionalOStream pcout;
};
template <int dim>
double coefficient(const Point<dim> &p)
{
return -2 * M_PI * M_PI * cos(M_PI * p[0]) * cos(M_PI * p[1]);
}
template <int dim>
Step5<dim>::Step5()
: mpi_communicator(MPI_COMM_WORLD),
triangulation(mpi_communicator,
typename Triangulation<dim>::MeshSmoothing(Triangulation<dim>::smoothing_on_refinement | Triangulation<dim>::smoothing_on_coarsening)),
dof_handler(triangulation),
fe(2),
pcout(std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
{}
template <int dim>
Step5<dim>::~Step5<dim>()
{
dof_handler.clear();
}
template <int dim>
void Step5<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
locally_owned_dofs = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
locally_relevant_solution.reinit(locally_owned_dofs, locally_relevant_dofs, mpi_communicator);
evaluation_point.reinit(locally_owned_dofs, locally_relevant_dofs, mpi_communicator);
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
hanging_node_constraints.clear();
hanging_node_constraints.reinit(locally_relevant_dofs);
VectorTools::interpolate_boundary_values(dof_handler,
1,
BoundaryValues<dim>(1),
hanging_node_constraints);
hanging_node_constraints.close();
pcout << " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
DynamicSparsityPattern dsp(locally_relevant_dofs);
DoFTools::make_sparsity_pattern(dof_handler, dsp, hanging_node_constraints, false);
SparsityTools::distribute_sparsity_pattern(dsp,
dof_handler.n_locally_owned_dofs_per_processor(),
mpi_communicator,
locally_relevant_dofs);
system_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
}
template <int dim>
void Step5<dim>::assemble_rhs(double &return_value,
const double factor,
const double &val_u,
const Tensor<1, dim> &grad_u,
const Point<dim> p,
const double &fe_values_value,
const Tensor<1, dim> &fe_gradient_value,
const double JxW_value)
{
(void) val_u;
return_value += factor
* (coefficient(p)
* fe_values_value
+ grad_u
* fe_gradient_value
)
* JxW_value;
}
template <int dim>
void Step5<dim>::calculate_residual(const TrilinosWrappers::MPI::Vector &u, TrilinosWrappers::MPI::Vector &residual)
{
QGauss<dim> quadrature_formula(fe.degree + 1);
UpdateFlags cell_flags = update_values | update_gradients | update_quadrature_points | update_JxW_values;
FEValues<dim> fe_values(fe, quadrature_formula, cell_flags);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
residual.reinit(locally_owned_dofs, mpi_communicator);
std::vector<Tensor<1, dim>> grad_u(quadrature_formula.size());
std::vector<double> val_u(quadrature_formula.size());
TrilinosWrappers::MPI::Vector local_evaluation_point(locally_relevant_dofs, mpi_communicator);
local_evaluation_point = u;
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
const unsigned int n_q_points = quadrature_formula.size();
for (const auto &cell : dof_handler.active_cell_iterators())
{
if(cell->is_locally_owned())
{
cell_rhs = 0.;
fe_values.reinit(cell);
fe_values.get_function_values(local_evaluation_point, val_u);
fe_values.get_function_gradients(local_evaluation_point, grad_u);
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
assemble_rhs(cell_rhs(i),
-1,
val_u[q],
grad_u[q],
fe_values.quadrature_point(q),
fe_values.shape_value(i, q),
fe_values.shape_grad(i, q),
fe_values.JxW(q));
}
}
cell->get_dof_indices(local_dof_indices);
hanging_node_constraints.distribute_local_to_global(cell_rhs,
local_dof_indices,
residual);
}
}
hanging_node_constraints.distribute(residual);
residual.compress(VectorOperation::add);
}
template <int dim>
void Step5<dim>::solve_with_NOX()
{
const double offset = 1;
TrilinosWrappers::MPI::Vector offset_vector(locally_owned_dofs, mpi_communicator);
offset_vector = offset;
//return;
TrilinosWrappers::MPI::Vector locally_owned_solution(locally_owned_dofs, mpi_communicator);
locally_owned_solution = offset;
Epetra_Vector local_solution(Copy, locally_owned_solution.trilinos_vector(), 0);
Teuchos::RCP<Teuchos::ParameterList> nlParamsPtr =
Teuchos::rcp(new Teuchos::ParameterList);
Teuchos::ParameterList& nlParams = *nlParamsPtr.get();
// Set the nonlinear solver method
nlParams.set("Nonlinear Solver", "Line Search Based");
// Set the printing parameters in the "Printing" sublist
Teuchos::ParameterList& printParams = nlParams.sublist("Printing");
printParams.set("Output Precision", 3);
printParams.set("Output Processor", 0);
printParams.set("Output Information",
NOX::Utils::Warning);
// Create printing utilities
NOX::Utils utils(printParams);
// Sublist for line search
Teuchos::ParameterList& searchParams = nlParams.sublist("Line Search");
searchParams.set("Method", "Full Step");
// Sublist for direction
Teuchos::ParameterList& dirParams = nlParams.sublist("Direction");
dirParams.set("Method", "Newton");
Teuchos::ParameterList& newtonParams = dirParams.sublist("Newton");
newtonParams.set("Forcing Term Method", "Constant");
// Sublist for linear solver for the Newton method
Teuchos::ParameterList& lsParams = newtonParams.sublist("Linear Solver");
lsParams.set("Aztec Solver", "GMRES");
lsParams.set("Max Iterations", 800);
lsParams.set("Tolerance", 1e-4);
lsParams.set("Preconditioner", "None");
//lsParams.set("Preconditioner", "Ifpack");
lsParams.set("Max Age Of Prec", 5);
Teuchos::RCP<ProblemInterface<dim>> interface =
Teuchos::rcp(new ProblemInterface<dim>(std::bind(&Step5::calculate_residual,
this,
std::placeholders::_1,
std::placeholders::_2),
mpi_communicator,
locally_owned_dofs,
locally_relevant_dofs));
Teuchos::RCP<NOX::Epetra::MatrixFree> MF =
Teuchos::rcp(new NOX::Epetra::MatrixFree(printParams, interface, local_solution));
// Create the linear system
Teuchos::RCP<NOX::Epetra::Interface::Required> iReq = interface;
Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac = MF;
Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec =
interface;
Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> linSys =
Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams,
iReq, iJac, MF,
local_solution));
// Create the Group
Teuchos::RCP<NOX::Epetra::Group> grp =
Teuchos::rcp(new NOX::Epetra::Group(printParams, iReq, local_solution,
linSys));
// Create the convergence tests
Teuchos::RCP<NOX::StatusTest::NormF> absresid =
Teuchos::rcp(new NOX::StatusTest::NormF(1.0e-8));
Teuchos::RCP<NOX::StatusTest::NormF> relresid =
Teuchos::rcp(new NOX::StatusTest::NormF(*grp.get(), 1.0e-2));
Teuchos::RCP<NOX::StatusTest::NormUpdate> update =
Teuchos::rcp(new NOX::StatusTest::NormUpdate(1.0e-5));
Teuchos::RCP<NOX::StatusTest::NormWRMS> wrms =
Teuchos::rcp(new NOX::StatusTest::NormWRMS(1.0e-2, 1.0e-8));
Teuchos::RCP<NOX::StatusTest::Combo> converged =
Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::AND));
converged->addStatusTest(absresid);
converged->addStatusTest(relresid);
converged->addStatusTest(wrms);
converged->addStatusTest(update);
Teuchos::RCP<NOX::StatusTest::MaxIters> maxiters =
Teuchos::rcp(new NOX::StatusTest::MaxIters(20));
Teuchos::RCP<NOX::StatusTest::FiniteValue> fv =
Teuchos::rcp(new NOX::StatusTest::FiniteValue);
Teuchos::RCP<NOX::StatusTest::Combo> combo =
Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
combo->addStatusTest(fv);
combo->addStatusTest(converged);
combo->addStatusTest(maxiters);
Teuchos::RCP<Teuchos::ParameterList> finalParamsPtr = nlParamsPtr;
// Create the method
Teuchos::RCP<NOX::Solver::Generic> solver =
NOX::Solver::buildSolver(grp, combo, finalParamsPtr);
NOX::StatusTest::StatusType status;
{
status = solver->solve();
}
if (status == NOX::StatusTest::Converged)
utils.out() << "Test Passed!" << std::endl;
else {
utils.out() << "Nonlinear solver failed to converge!" << std::endl;
}
// Get the Epetra_Vector with the final solution from the solver
const NOX::Epetra::Group& finalGroup = dynamic_cast<const NOX::Epetra::Group&>(solver->getSolutionGroup());
const Epetra_Vector& finalSolution = (dynamic_cast<const NOX::Epetra::Vector&>(finalGroup.getX())).getEpetraVector();
// End Nonlinear Solver **************************************
// Output the parameter list
if (utils.isPrintType(NOX::Utils::Parameters)) {
utils.out() << std::endl << "Final Parameters" << std::endl
<< "****************" << std::endl;
solver->getList().print(utils.out());
utils.out() << std::endl;
}
locally_relevant_solution = 0;
Epetra_Vector epetra_F(View, locally_owned_solution.trilinos_vector(), 0);
epetra_F = finalSolution;
locally_owned_solution -= offset_vector;
//hanging_node_constraints.distribute(locally_owned_solution);
locally_relevant_solution = locally_owned_solution;
//hanging_node_constraints.distribute(locally_relevant_solution);
}
template <int dim>
void Step5<dim>::run()
{
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(2);
setup_system();
solve_with_NOX();
}
int main(int argc, char *argv[])
{
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
Step5<2> laplace_problem_2d;
laplace_problem_2d.run();
return 0;
}
