From 9543118d0d3dc8666c4e910b308caba2ab724251 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20Dutheillet-Lamonth=C3=A9zie?= Date: Wed, 23 Aug 2023 16:15:57 +0200 Subject: [PATCH] #132: start Laplace2D header for Tpetra --- .../tpetra/NOX_Operators/Laplace2D_Tpetra.hpp | 224 ++++++++++++++++++ 1 file changed, 224 insertions(+) create mode 100644 packages/nox/test/tpetra/NOX_Operators/Laplace2D_Tpetra.hpp diff --git a/packages/nox/test/tpetra/NOX_Operators/Laplace2D_Tpetra.hpp b/packages/nox/test/tpetra/NOX_Operators/Laplace2D_Tpetra.hpp new file mode 100644 index 000000000000..475782e88198 --- /dev/null +++ b/packages/nox/test/tpetra/NOX_Operators/Laplace2D_Tpetra.hpp @@ -0,0 +1,224 @@ +//@HEADER +// ************************************************************************ +// +// NOX: An Object-Oriented Nonlinear Solver Package +// Copyright (2002) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger Pawlowski (rppawlo@sandia.gov) or +// Eric Phipps (etphipp@sandia.gov), Sandia National Laboratories. +// ************************************************************************ +// CVS Information +// $Source$ +// $Author$ +// $Date$ +// $Revision$ +// ************************************************************************ +//@HEADER + +// NOX Test Lacplace 2D Problem +// ---------------------------- +// Simple nonlinear PDE problem. +// This file provides the problem +// +// -\Delta u + \lambda e^u = 0 in \Omega = (0,1) \times (0,1) +// u = 0 on \partial \Omega +// +// for use as a driver for various nox-epetra tests + +#include "Tpetra_Map.hpp" +#include "Tpetra_Vector.hpp" +#include "Tpetra_CrsMatrix.hpp" +#include "NOX_TpetraTypedefs.hpp" +#include "NOX.H" +// #include "NOX_Epetra_Interface_Required.H" +// #include "NOX_Epetra_Interface_Jacobian.H" +// #include "NOX_Epetra_Interface_Preconditioner.H" +// #include "NOX_Epetra_LinearSystem_AztecOO.H" + +// Typedefs + +// typedef Tpetra::Vector<>::scalar_type Scalar; +// typedef Tpetra::Vector<>::local_ordinal_type LO; +// typedef Tpetra::Vector<>::global_ordinal_type GO; +// typedef Tpetra::Vector<>::node_type Node; +typedef NOX::Scalar Scalar; +typedef NOX::LocalOrdinal LO; +typedef NOX::GlobalOrdinal GO; +typedef NOX::NodeType Node; + +typedef Tpetra::Map Tpetra_Map; +typedef Tpetra::Vector Tpetra_Vector; +typedef Tpetra::MultiVector Tpetra_MultiVector; +typedef Tpetra::CrsMatrix Tpetra_CsrMatrix; +typedef Tpetra::Operator Tpetra_Operator +// typedef Thyra::VectorSpaceBase TVSB; +// typedef Thyra::VectorBase TVB; +// typedef Thyra::MultiVectorBase TMVB; +// typedef Thyra::LinearOpBase TLOB; +// typedef Thyra::TpetraOperatorVectorExtraction TOVE; +// typedef NOX::Thyra::Vector NTV; +// typedef NOX::Thyra::MultiVector NTMV; +// typedef typename TMV::mag_type mag_type; + +namespace Laplace2D { + +// this is required to know the number of lower, upper, left and right +// node for each node of the Cartesian grid (composed by nx \timex ny +// elements) + +void +get_myNeighbours( const int i, const int nx, const int ny, + int & left, int & right, + int & lower, int & upper); + + +// This function creates a CrsMatrix, whose elements corresponds +// to the discretization of a Laplacian over a Cartesian grid, +// with nx grid point along the x-axis and and ny grid points +// along the y-axis. For the sake of simplicity, I suppose that +// all the nodes in the matrix are internal nodes (Dirichlet +// boundary nodes are supposed to have been already condensated) + +Tpetra_CsrMatrix * +CreateLaplacian( const int nx, const int ny, const Teuchos::RCP >& comm); + +// ========================================================================== +// This class contians the main definition of the nonlinear problem at +// hand. A method is provided to compute F(x) for a given x, and another +// method to update the entries of the Jacobian matrix, for a given x. +// As the Jacobian matrix J can be written as +// J = L + diag(lambda*exp(x[i])), +// where L corresponds to the discretization of a Laplacian, and diag +// is a diagonal matrix with lambda*exp(x[i]). Basically, to update +// the jacobian we simply update the diagonal entries. Similarly, to compute +// F(x), we reset J to be equal to L, then we multiply it by the +// (distributed) vector x, then we add the diagonal contribution +// ========================================================================== + +} // namespace Laplace2D + + +class PDEProblem +{ + +public: + + // constructor. Requires the number of nodes along the x-axis + // and y-axis, the value of lambda, and the Epetra_Communicator + // (to define a Map, which is a linear map in this case) + PDEProblem(const int nx, const int ny, const double lambda, + const Teuchos::RCP >& comm); + + // destructor + ~PDEProblem(); + + // compute F(x) + void ComputeF(const TV & x, TV & f); + + // update the Jacobian matrix for a given x + void UpdateJacobian(const TV & x); + + // returns a pointer to the internally stored matrix + Tpetra_CsrMatrix * GetMatrix() + { + return Matrix_; + } + +private: + + int nx_, ny_; + double hx_, hy_; + Tpetra_CsrMatrix * Matrix_; + double lambda_; + +}; /* class PDEProblem */ + +// ========================================================================== +// This is the main NOX class for this example. Here we define +// the interface between the nonlinear problem at hand, and NOX. +// The constructor accepts a PDEProblem object. Using a pointer +// to this object, we can update the Jacobian and compute F(x), +// using the definition of our problem. This interface is bit +// crude: For instance, no PrecMatrix nor Preconditioner is specified. +// ========================================================================== + +class SimpleProblemInterface : public NOX::Epetra::Interface::Required , // TODO: no tpetra alternaative create it + public NOX::Epetra::Interface::Jacobian , // TODO: no tpetra alternaative create it + public NOX::Epetra::Interface::Preconditioner // TODO: no tpetra alternaative create it +{ + +public: + + //! Constructor + SimpleProblemInterface( PDEProblem * Problem ) : + Problem_(Problem) {}; + + //! Destructor + ~SimpleProblemInterface() + { + }; + + bool computeF(const Tpetra_Vector & x, Tpetra_Vector & f, + NOX::Epetra::Interface::Required::FillType F ) + { + Problem_->ComputeF(x,f); + return true; + }; + + bool computeJacobian(const Tpetra_Vector & x, Tpetra_Operator & Jac) + { + + Problem_->UpdateJacobian(x); + return true; + } + + bool computePreconditioner(const Tpetra_Vector & x, Tpetra_Operator & Op, Teuchos::ParameterList *) + { + + Problem_->UpdateJacobian(x); + return true; + } + + bool computePrecMatrix(const Tpetra_Vector & x, Tpetra_RowMatrix & M) + { + std::cout << "*ERR* SimpleProblem::preconditionVector()\n"; + std::cout << "*ERR* don't use explicit preconditioning" << std::endl; + throw 1; + } + +private: + + PDEProblem * Problem_; + +}; /* class SimpleProblemInterface */ +