diff --git a/packages/nox/test/tpetra/CMakeLists.txt b/packages/nox/test/tpetra/CMakeLists.txt index 5334a5313c19..ef652a2ddaab 100644 --- a/packages/nox/test/tpetra/CMakeLists.txt +++ b/packages/nox/test/tpetra/CMakeLists.txt @@ -68,3 +68,6 @@ IF(NOX_ENABLE_ABSTRACT_IMPLEMENTATION_THYRA AND ENDIF() ENDIF() + +ADD_SUBDIRECTORY(LOCA_TestProblems) +ADD_SUBDIRECTORY(LOCA_UnitTests) diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp b/packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp new file mode 100644 index 000000000000..530776a281dd --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp @@ -0,0 +1,89 @@ +// $Id$ +// $Source$ + +//@HEADER +// ************************************************************************ +// +// LOCA: Library of Continuation Algorithms Package +// Copyright (2005) 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 + +#include "Basis.hpp" + +// Constructor +Basis::Basis() : + phi(2), + dphide(2), + uu(0.0), + xx(0.0), + duu(0.0), + eta(0.0), + wt(0.0), + dx(0.0) +{ } + +// Calculates a linear 1D basis +void Basis::getBasis(int gp, std::vector x, std::vector u) { + int N = 2; + if (gp==0) {eta=-1.0/std::sqrt(3.0); wt=1.0;} + if (gp==1) {eta=1.0/std::sqrt(3.0); wt=1.0;} + + // Calculate basis function and derivatives at nodel pts + phi[0]=(1.0-eta)/2.0; + phi[1]=(1.0+eta)/2.0; + dphide[0]=-0.5; + dphide[1]=0.5; + + // Caculate basis function and derivative at GP. + dx=0.5*(x[1]-x[0]); + xx=0.0; + uu=0.0; + duu=0.0; + for (int i=0; i < N; i++) { + xx += x[i] * phi[i]; + uu += u[i] * phi[i]; + duu += u[i] * dphide[i]; + } + + return; +} diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/Basis.hpp b/packages/nox/test/tpetra/LOCA_TestProblems/Basis.hpp new file mode 100644 index 000000000000..1f76a0a8f5cf --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/Basis.hpp @@ -0,0 +1,81 @@ +// $Id$ +// $Source$ + +//@HEADER +// ************************************************************************ +// +// LOCA: Library of Continuation Algorithms Package +// Copyright (2005) 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 + +// 1D Linear basis function for finite element method + +#ifndef _NOX_EXAMPLE_TPETRA_LINEAR_BASIS_H +#define _NOX_EXAMPLE_TPETRA_LINEAR_BASIS_H + +#include +#include +#include + +class Basis +{ + public: + // Constructor + Basis(); + + // Calculates the values of u and x at the specified gauss point + void getBasis(int gp, std::vector x, std::vector u); + + private: + // Private to prohibit copying + Basis(const Basis&); + Basis& operator=(const Basis&); + + public: + // Variables that are calculated at the gauss point + std::vector phi, dphide; + double uu, xx, duu, eta, wt; + double dx; +}; + +#endif diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt new file mode 100644 index 000000000000..62035204ecee --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt @@ -0,0 +1,39 @@ + +SET(HEADERS "") +SET(SOURCES "") + +TRIBITS_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}) + +APPEND_SET(HEADERS + Basis.hpp + FiniteElementProblem.hpp +# Problem_Interface.hpp +# Tcubed_FiniteElementProblem.hpp + Pitchfork_FiniteElementProblem.hpp # --> WIP, TESTING +# NormConstraint.hpp +# LOCALinearConstraint.hpp --> WIP +) + +APPEND_SET(SOURCES + Basis.cpp +# Problem_Interface.cpp +# Tcubed_FiniteElementProblem.cpp + Pitchfork_FiniteElementProblem.cpp # --> WIP, TESTING +# NormConstraint.cpp +# LOCALinearConstraint.cpp --> WIP +) + +# Add this below: +# IF(NOX_ENABLE_ABSTRACT_IMPLEMENTATION_EPETRA AND NOX_ENABLE_LOCA) + +IF(NOX_ENABLE_LOCA) + + TRIBITS_ADD_LIBRARY( + locatpetratestproblems + HEADERS ${HEADERS} + SOURCES ${SOURCES} + TESTONLY + DEPLIBS loca locatpetra + ) + +ENDIF() diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp new file mode 100644 index 000000000000..a2ac87da976b --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp @@ -0,0 +1,310 @@ +//@HEADER +// ************************************************************************ +// +// LOCA: Library of Continuation Algorithms Package +// Copyright (2005) 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 + +#include + +#include "Basis.hpp" + +#include "FiniteElementProblem.hpp" + +// Constructor - creates the Tpetra objects (maps and vectors) +template +FiniteElementProblem::FiniteElementProblem(int numGlobalElements, + Teuchos::RCP>& comm, + Scalar s) + : comm(&comm), numGlobalElements(numGlobalElements), scale(s) { + // Commonly used variables + int i; + myRank = comm->getRank(); // Process ID + numProc = comm->getSize(); // Total number of processes + + // Construct a Source Map that puts approximately the same + // Number of equations on each processor in uniform global ordering + standardMap = new map_type(numGlobalElements, 0, *comm); + + // Get the number of elements owned by this processor + numMyElements = standardMap->getLocalNumElements(); + + // Construct an overlaped map for the finite element fill ************* + // For single processor jobs, the overlap and standard map are the same + if (numProc == 1) { + overlapMap = new map_type(*standardMap); + } else { + int overlapNumMyElements; + int overlapMinMyGID; + overlapNumMyElements = numMyElements + 2; + if ((myRank == 0) || (myRank == numProc - 1)) + overlapNumMyElements--; + + if (myRank == 0) { + overlapMinMyGID = standardMap->getMinGlobalIndex(); + } else { + overlapMinMyGID = standardMap->getMinGlobalIndex() - 1; + } + + int* overlapMyGlobalElements = new int[overlapNumMyElements]; + + for (i = 0; i < overlapNumMyElements; i++) { + overlapMyGlobalElements[i] = overlapMinMyGID + i; + } + + overlapMap = new tmap_t(-1, overlapNumMyElements, overlapMyGlobalElements, 0, *comm); + + } // End Overlap map construction ************************************* + + // Construct Linear Objects + importer = new timport_t(*overlapMap, *standardMap); + initialSolution = new tvector_t(*standardMap); + AA = new tcrsgraph_t(Copy, *standardMap, 5); + + // Allocate the memory for a matrix dynamically (i.e. the graph is dynamic). + generateGraph(*AA); + + // Create a second matrix using graph of first matrix - this creates a + // static graph so we can refill the new matirx after FillComplete() + // is called. + A = new tcrsmatrix_t(Copy, *AA); + A->FillComplete(); + + // Set default bifurcation values + factor = 0.1; + leftBC = 1.0; + rightBC = 1.0; +} + +// Matrix and Residual Fills +template +bool FiniteElementProblem::evaluate(FillType f, const tvector_t* soln, tvector_t* tmp_rhs, + trowmatrix_t* tmp_matrix) { + flag = f; + + // Set the incoming linear objects + if (flag == F_ONLY) { + rhs = tmp_rhs; + } else if (flag == MATRIX_ONLY) { + A = dynamic_cast(tmp_matrix); + } else if (flag == ALL) { + rhs = tmp_rhs; + A = dynamic_cast(tmp_matrix); + } else { + std::cout << "ERROR: FiniteElementProblem::fillMatrix() - FillType flag is broken" << std::endl; + throw; + } + + // Create the overlapped solution and position vectors + vector_type u(*overlapMap); + vector_type x(*overlapMap); + + // Export Solution to Overlap vector + u.doImport(*soln, *importer, Tpetra::INSERT); + + // Declare required variables + int i, j, ierr; + int overlapNumMyElements = overlapMap->numMyElements(); + + int overlapMinMyGID; + if (myRank == 0) + overlapMinMyGID = standardMap->getMinGlobalIndex(); + else + overlapMinMyGID = standardMap->getMinGlobalIndex() - 1; + + int row, column; + ST jac; + ST xx[2]; + ST uu[2]; + Basis basis; + + // Create the nodal coordinates + ST Length = 1.0; + ST dx = Length / ((ST)numGlobalElements - 1); + for (i = 0; i < overlapNumMyElements; i++) { + x[i] = dx * ((ST)overlapMinMyGID + i); + } + + // Zero out the objects that will be filled + if ((flag == MATRIX_ONLY) || (flag == ALL)) + i = A->PutScalar(0.0); + if ((flag == F_ONLY) || (flag == ALL)) + i = rhs->PutScalar(0.0); + + // Loop Over # of Finite Elements on Processor + for (int ne = 0; ne < overlapNumMyElements - 1; ne++) { + // Loop Over Gauss Points + for (int gp = 0; gp < 2; gp++) { + // Get the solution and coordinates at the nodes + xx[0] = x[ne]; + xx[1] = x[ne + 1]; + uu[0] = u[ne]; + uu[1] = u[ne + 1]; + // Calculate the basis function at the gauss point + basis.getBasis(gp, xx, uu); + + // Loop over Nodes in Element + for (i = 0; i < 2; i++) { + row = overlapMap->getGlobalElement(ne + i); + // printf("Proc=%d GlobalRow=%d LocalRow=%d Owned=%d\n", + // myRank, row, ne+i,standardMap.MyGID(row)); + if (standardMap->MyGID(row)) { + if ((flag == F_ONLY) || (flag == ALL)) { + (*rhs)[standardMap->LID(overlapMap->GID(ne + i))] += + +basis.wt * basis.dx * + ((-1.0 / (basis.dx * basis.dx)) * basis.duu * basis.dphide[i] + + factor * basis.uu * basis.uu * basis.uu * basis.phi[i]); + } + } + // Loop over Trial Functions + if ((flag == MATRIX_ONLY) || (flag == ALL)) { + for (j = 0; j < 2; j++) { + if (standardMap->getLocalElement(row)) { + column = overlapMap->GID(ne + j); + jac = basis.wt * basis.dx * + ((-1.0 / (basis.dx * basis.dx)) * basis.dphide[j] * basis.dphide[i] + + 3.0 * factor * basis.uu * basis.uu * basis.phi[j] * basis.phi[i]); + ierr = A->sumIntoGlobalValues(row, 1, &jac, &column); + } + } + } + } + } + } + + // Insert Boundary Conditions and modify Jacobian and function (F) + // U(0)=1 + if (myRank == 0) { + if ((flag == F_ONLY) || (flag == ALL)) + (*rhs)[0] = (*soln)[0] - leftBC; + if ((flag == MATRIX_ONLY) || (flag == ALL)) { + column = 0; + jac = 1.0; + A->replaceGlobalValues(0, 1, &jac, &column); + column = 1; + jac = 0.0; + A->replaceGlobalValues(0, 1, &jac, &column); + } + } + + if (standardMap->LID(standardMap->getMaxAllGlobalIndex()) >= 0) { + int lastDof = standardMap->getLocalElement(standardMap->getMaxAllGlobalIndex()); + if ((flag == F_ONLY) || (flag == ALL)) + (*rhs)[lastDof] = (*soln)[lastDof] - rightBC; + if ((flag == MATRIX_ONLY) || (flag == ALL)) { + int row = standardMap->getMaxAllGlobalIndex(); + column = row; + jac = 1.0; + A->replaceGlobalValues(row, 1, &jac, &column); + column = row - 1; + jac = 0.0; + A->replaceGlobalValues(row, 1, &jac, &column); + } + } + // Sync up processors to be safe + comm->barrier(); + + A->fillComplete(); + + return true; +} + +template +Teuchos::RCP FiniteElementProblem::getSolution() { + return initialSolution; +} + +template +Teuchos::RCP FiniteElementProblem::getJacobian() { + return A; +} + +template +bool FiniteElementProblem::setParameter(std::string label, ST value) { + if (label == "Nonlinear Factor") + factor = value; + else if (label == "Left BC") + leftBC = value; + else if (label == "Right BC") + rightBC = value * scale; + else if (label == "Homotopy Continuation Parameter") { + // do nothing for now + } else { + std::cout << "ERROR: FiniteElementProblem::setParameter() - label is invalid " + << "for this problem!" << std::endl; + exit(-1); + } + return true; +} + +template +tcrsgraph_t& FiniteElementProblem::generateGraph(tcrsgraph_t& AAA) { + // Declare required variables + int i, j; + int row, column; + int overlapNumMyElements = overlapMap->getLocalNumElements(); + int overlapMinMyGID; + if (myRank == 0) + overlapMinMyGID = standardMap->getMinGlobalIndex(); + else + overlapMinMyGID = standardMap->getMinGlobalIndex() - 1; + + // Loop Over # of Finite Elements on Processor + for (int ne = 0; ne < overlapNumMyElements - 1; ne++) { + // Loop over Nodes in Element + for (i = 0; i < 2; i++) { + row = overlapMap->GID(ne + i); + + // Loop over Trial Functions + for (j = 0; j < 2; j++) { + // If this row is owned by current processor, add the index + if (standardMap->MyGID(row)) { + column = overlapMap->GID(ne + j); + AAA.insertGlobalIndices(row, 1, &column); + } + } + } + } + AAA.fillComplete(); + return AAA; +} diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.hpp b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.hpp new file mode 100644 index 000000000000..33c943f8389f --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.hpp @@ -0,0 +1,109 @@ +// $Id$ +// $Source$ + +//@HEADER +// ************************************************************************ +// +// LOCA: Library of Continuation Algorithms Package +// Copyright (2005) 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 + +// Finite Element Problem Class + +#ifndef _FINITEELEMENTPROBLEM_H +#define _FINITEELEMENTPROBLEM_H + +#include + +#include +#include +#include +#include +#include +#include + +// Flag to tell the evaluate routine what objects to fill +enum FillType {F_ONLY, MATRIX_ONLY, ALL}; + +// Finite Element Problem Class +template +class FiniteElementProblem { + + using LO = typename Tpetra::Vector<>::local_ordinal_type; + using GO = typename Tpetra::Vector<>::global_ordinal_type; + using NO = typename Tpetra::Map<>::node_type; + +public: + + using scalar_type = Scalar; + using vector_type = typename Tpetra::Vector; + using row_matrix_type = typename Tpetra::RowMatrix; + using crs_matrix_type = typename Tpetra::CrsMatrix; + using import_type = typename Tpetra::Import; + using csr_graph_type = typename Tpetra::CrsGraph; + using map_type = typename Tpetra::Map; + + // Constructor + FiniteElementProblem() {} + + // Evaluates the function (F) and/or the Jacobian using the solution + // values in solnVector. + virtual bool evaluate(FillType f, const vector_type *solnVector, + vector_type *rhsVector, + row_matrix_type *matrix, + scalar_type jac_coeff = 1.0, + scalar_type mass_coeff = 0.0) = 0; + + // Return a reference to the Tpetra_Vector with the initial guess + // that is generated by the FiniteElementProblem class. + virtual vector_type& getSolution() = 0; + + // Return a reference to the Tpetra_Vector with the Jacobian + // that is generated by the FiniteElementProblem class. + virtual crs_matrix_type& getJacobian() = 0; + + // Set a bifurcation parameter in the application physics + virtual bool setParameter(std::string label, scalar_type value) = 0; + +}; +#endif diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/LOCALinearConstraint.cpp b/packages/nox/test/tpetra/LOCA_TestProblems/LOCALinearConstraint.cpp new file mode 100644 index 000000000000..b706a5f71952 --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/LOCALinearConstraint.cpp @@ -0,0 +1,230 @@ +// $Id$ +// $Source$ + +//@HEADER +// ************************************************************************ +// +// LOCA: Library of Continuation Algorithms Package +// Copyright (2005) 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 + +#include "LOCALinearConstraint.H" +#include "LOCA_Parameter_Vector.H" + +LinearConstraint::LinearConstraint(int NumConstraints, + const LOCA::ParameterVector& pVec, + NOX::Abstract::Vector& cloneVec) : + m(NumConstraints), + constraints(m,1), + isValidConstraints(false), + x(), + dgdx(), + isZeroDgDx(false), + p(pVec), + pvec(Teuchos::View, p.getDoubleArrayPointer(), p.length(), p.length(), 1), + dgdp(NumConstraints, pVec.length()) +{ + constraints.putScalar(0.0); + x = cloneVec.createMultiVector(1); + dgdx = cloneVec.createMultiVector(m); + dgdp.putScalar(0.0); +} + +LinearConstraint::LinearConstraint(const LinearConstraint& source, + NOX::CopyType type) : + m(source.m), + constraints(source.constraints), + isValidConstraints(false), + x(source.x->clone(type)), + dgdx(source.dgdx->clone(type)), + isZeroDgDx(source.isZeroDgDx), + p(source.p), + pvec(source.pvec), + dgdp(source.dgdp) +{ + if (source.isValidConstraints && type == NOX::DeepCopy) + isValidConstraints = true; +} + +LinearConstraint::~LinearConstraint() +{ +} + +void +LinearConstraint::copy(const LOCA::MultiContinuation::ConstraintInterface& src) +{ + const LinearConstraint& source = dynamic_cast(src); + + if (this != &source) { + m = source.m; + constraints = source.constraints; + isValidConstraints = source.isValidConstraints; + *x = *source.x; + *dgdx = *source.dgdx; + isZeroDgDx = source.isZeroDgDx; + p = source.p; + pvec.assign(source.pvec); + dgdp = source.dgdp; + } +} + +Teuchos::RCP +LinearConstraint::clone(NOX::CopyType type) const +{ + return Teuchos::rcp(new LinearConstraint(*this, type)); +} + +int +LinearConstraint::numConstraints() const +{ + return m; +} + +void +LinearConstraint::setX(const NOX::Abstract::Vector& y) +{ + (*x)[0] = y; + isValidConstraints = false; +} + +void +LinearConstraint::setParam(int paramID, double val) +{ + p[paramID] = val; +} + +void +LinearConstraint::setParams( + const std::vector& paramIDs, + const NOX::Abstract::MultiVector::DenseMatrix& vals) +{ + for (unsigned int i=0; imultiply(1.0, *dgdx, constraints); + if (pvec.numRows() > 0) + constraints.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, dgdp, + pvec, 1.0); + isValidConstraints = true; + + return NOX::Abstract::Group::Ok; +} + +NOX::Abstract::Group::ReturnType +LinearConstraint::computeDX() +{ + return NOX::Abstract::Group::Ok; +} + +NOX::Abstract::Group::ReturnType +LinearConstraint::computeDP(const std::vector& paramIDs, + NOX::Abstract::MultiVector::DenseMatrix& dp, + bool isValidG) +{ + if (!isValidG) { + for (int i=0; i + clone(NOX::CopyType type = NOX::DeepCopy) const; + + //! Return number of constraints + virtual int numConstraints() const; + + //! Set the solution vector to y. + virtual void setX(const NOX::Abstract::Vector& y); + + //! Sets parameter indexed by paramID + virtual void setParam(int paramID, double val); + + //! Sets parameters indexed by paramIDs + virtual void setParams(const std::vector& paramIDs, + const NOX::Abstract::MultiVector::DenseMatrix& vals); + + //! Compute continuation constraint equations + virtual NOX::Abstract::Group::ReturnType + computeConstraints(); + + //! Compute derivative of constraints w.r.t. solution vector x + virtual NOX::Abstract::Group::ReturnType + computeDX(); + + //! Compute derivative of constraints w.r.t. supplied parameters. + /*! + * The first column of \c dgdp should be filled with the constraint + * residuals \f$g\f$ if \c isValidG is \c false. If \c isValidG is + * \c true, then the \c dgdp contains \f$g\f$ on input. + */ + virtual NOX::Abstract::Group::ReturnType + computeDP(const std::vector& paramIDs, + NOX::Abstract::MultiVector::DenseMatrix& dp, + bool isValidG); + + //! Return \c true if constraint residuals are valid + virtual bool isConstraints() const; + + //! Return \c true if derivative of constraint w.r.t. x is valid + virtual bool isDX() const; + + //! Return constraint residuals + virtual const NOX::Abstract::MultiVector::DenseMatrix& + getConstraints() const; + + //! Return solution component of constraint derivatives + virtual const NOX::Abstract::MultiVector* + getDX() const; + + /*! + * \brief Return \c true if solution component of constraint + * derivatives is zero + */ + virtual bool isDXZero() const; + + //@} + + //! Set solution component of constraint derivative + void setDgDx(const NOX::Abstract::MultiVector& A); + + //! Set parameter component of constraint derivative + void setDgDp(const NOX::Abstract::MultiVector::DenseMatrix& A); + + /*! + * \brief Set the flag determinging whether isDXZero() + * returns true or false. This is for testing bordered solvers when B + * is zero. + */ + void setIsZeroDX(bool flag); + +private: + + //! Prohibit generation and use of operator=() + LinearConstraint& operator=(const LinearConstraint& source); + +protected: + + //! Number of constraints + int m; + + //! Constraint values + NOX::Abstract::MultiVector::DenseMatrix constraints; + + //! Flag indicating whether constraints are valid + bool isValidConstraints; + + //! Solution vector + Teuchos::RCP x; + + //! Solution component of constraint derivative + Teuchos::RCP dgdx; + + //! Flag indicating whether dgdx is zero + bool isZeroDgDx; + + //! Parameter vector + LOCA::ParameterVector p; + + //! Parameter vector as a DenseMatrix (view of p) + NOX::Abstract::MultiVector::DenseMatrix pvec; + + //! Parameter component of constraint derivative + NOX::Abstract::MultiVector::DenseMatrix dgdp; + +}; + +#endif diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/Pitchfork_FiniteElementProblem.cpp b/packages/nox/test/tpetra/LOCA_TestProblems/Pitchfork_FiniteElementProblem.cpp new file mode 100644 index 000000000000..29bccda3735c --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/Pitchfork_FiniteElementProblem.cpp @@ -0,0 +1,347 @@ +// $Id$ +// $Source$ + +//@HEADER +// ************************************************************************ +// +// LOCA: Library of Continuation Algorithms Package +// Copyright (2005) 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 + +#include "Pitchfork_FiniteElementProblem.hpp" + +#include + +#include "Basis.H" +#include "NOX_Common.H" + +// Constructor - creates the Tpetra objects (maps and vectors) +template +Pitchfork_FiniteElementProblem::Pitchfork_FiniteElementProblem( + int numGlobalElements, Teuchos::RCP>& comm) + : flag(F_ONLY), + standardMap(NULL), + overlapMap(NULL), + importer(NULL), + initialSolution(NULL), + rhs(NULL), + AA(NULL), + A(NULL), + comm(&comm), + rank(0), + numProcs(0), + numLocalElements(0), + numGlobalElements(numGlobalElements), + lambda(0.0), + alpha(0.0), + beta(0.0) { + // Commonly used variables + int i; + rank = comm->getRank(); // Process ID + numProcs = comm->NumProc(); // Total number of processes + + // Construct a Source Map that puts approximately the same + // Number of equations on each processor in uniform global ordering + standardMap = new map_type(numGlobalElements, 0, *comm); + + // Get the number of elements owned by this processor + numLocalElements = standardMap->getLocalNumElements(); + + // Construct an overlaped map for the finite element fill ************* + // For single processor jobs, the overlap and standard map are the same + if (NumProc == 1) { + overlapMap = new map_type(*standardMap); + } else { + int overlapNumLocalElements; + int overlapMinMyGID; + + overlapNumLocalElements = numLocalElements + 2; + if ((myRank == 0) || (myRank == NumProc - 1)) + overlapNumLocalElements--; + + if (myRank == 0) { + overlapMinMyGID = standardMap->MinMyGID(); + } + else { + overlapMinMyGID = standardMap->MinMyGID() - 1; + } + + int* overlapMyGlobalElements = new int[overlapNumLocalElements]; + + for (i = 0; i < overlapNumLocalElements; i++) { + overlapMyGlobalElements[i] = overlapMinMyGID + i; + } + + overlapMap = new map_type(-1, overlapNumLocalElements, OverlapMyGlobalElements, 0, *comm); + + delete[] overlapMyGlobalElements; + + } // End Overlap map construction ************************************* + + // Construct Linear Objects + importer = new import_type(*overlapMap, *standardMap); + initialSolution = new vector_type(*standardMap); + AA = new csr_graph_type(Copy, *standardMap, 5); + + // Allocate the memory for a matrix dynamically (i.e. the graph is dynamic). + generateGraph(*AA); + + // Create a second matrix using graph of first matrix - this creates a + // static graph so we can refill the new matirx after FillComplete() + // is called. + A = new csr_matrix_type(Copy, *AA); + A->fillComplete(); + + // Set default bifurcation values + lambda = -2.25; + alpha = 1.0; + beta = 0.0; +} + +// Destructor +template +Pitchfork_FiniteElementProblem::~Pitchfork_FiniteElementProblem() { + delete AA; + delete A; + delete initialSolution; + delete importer; + delete overlapMap; + delete standardMap; +} + +// Matrix and Residual Fills +template +bool Pitchfork_FiniteElementProblem::evaluate(FillType f, const Tpetra_Vector* soln, + Tpetra_Vector* tmp_rhs, Epetra_RowMatrix* tmp_matrix, + double jac_coeff, double mass_coeff) { + flag = f; + + // Set the incoming linear objects + if (flag == F_ONLY) { + rhs = tmp_rhs; + } else if (flag == MATRIX_ONLY) { + A = dynamic_cast(tmp_matrix); + assert(A != NULL); + } else if (flag == ALL) { + rhs = tmp_rhs; + A = dynamic_cast(tmp_matrix); + assert(A != NULL); + } else { + std::cout << "ERROR: Pitchfork_FiniteElementProblem::fillMatrix() - FillType flag is broken" + << std::endl; + throw; + } + + // Create the overlapped solution and position vectors + Tpetra_Vector u(*overlapMap); + Tpetra_Vector x(*overlapMap); + + // Export Solution to Overlap vector + u.Import(*soln, *Importer, Insert); + + // Declare required variables + int i, j, ierr; + int overlapNumLocalElements = overlapMap->NumMyElements(); + + int overlapMinMyGID; + if (myRank == 0) + overlapMinMyGID = standardMap->MinMyGID(); + else + overlapMinMyGID = standardMap->MinMyGID() - 1; + + int row, column; + double jac; + double xx[2]; + double uu[2]; + Basis basis; + + // Create the nodal coordinates + double Length = 2.0; + double dx = Length / ((double)NumGlobalElements - 1); + for (i = 0; i < overlapNumLocalElements; i++) { + x[i] = -1.0 + dx * ((double)overlapMinMyGID + i); + } + + // Zero out the objects that will be filled + if ((flag == MATRIX_ONLY) || (flag == ALL)) { + i = A->PutScalar(0.0); + assert(i == 0); + } + if ((flag == F_ONLY) || (flag == ALL)) { + i = rhs->PutScalar(0.0); + assert(i == 0); + } + + // Loop Over # of Finite Elements on Processor + for (int ne = 0; ne < overlapNumLocalElements - 1; ne++) { + // Loop Over Gauss Points + for (int gp = 0; gp < 2; gp++) { + // Get the solution and coordinates at the nodes + xx[0] = x[ne]; + xx[1] = x[ne + 1]; + uu[0] = u[ne]; + uu[1] = u[ne + 1]; + // Calculate the basis function at the gauss point + basis.getBasis(gp, xx, uu); + + // Loop over Nodes in Element + for (i = 0; i < 2; i++) { + row = overlapMap->GID(ne + i); + // printf("Proc=%d GlobalRow=%d LocalRow=%d Owned=%d\n", + // myRank, row, ne+i,standardMap.MyGID(row)); + if (standardMap->MyGID(row)) { + if ((flag == F_ONLY) || (flag == ALL)) { + (*rhs)[standardMap->LID(overlapMap->GID(ne + i))] += + +basis.wt * basis.dx * + ((-1.0 / (basis.dx * basis.dx)) * basis.duu * basis.dphide[i] - + source_term(basis.uu) * basis.phi[i]); + } + } + // Loop over Trial Functions + if ((flag == MATRIX_ONLY) || (flag == ALL)) { + for (j = 0; j < 2; j++) { + if (standardMap->MyGID(row)) { + column = overlapMap->GID(ne + j); + jac = jac_coeff * basis.wt * basis.dx * + ((-1.0 / (basis.dx * basis.dx)) * basis.dphide[j] * basis.dphide[i] - + source_deriv(basis.uu) * basis.phi[j] * basis.phi[i]) + + mass_coeff * basis.wt * basis.dx * basis.phi[j] * basis.phi[i]; + ierr = A->SumIntoGlobalValues(row, 1, &jac, &column); + TEUCHOS_ASSERT(ierr == 0); + } + } + } + } + } + } + + // Insert Boundary Conditions and modify Jacobian and function (F) + // U(-1)=beta + if (myRank == 0) { + if ((flag == F_ONLY) || (flag == ALL)) + (*rhs)[0] = (*soln)[0] - beta; + if ((flag == MATRIX_ONLY) || (flag == ALL)) { + column = 0; + jac = 1.0 * jac_coeff; + A->ReplaceGlobalValues(0, 1, &jac, &column); + column = 1; + jac = 0.0 * jac_coeff; + A->ReplaceGlobalValues(0, 1, &jac, &column); + } + } + + // U(1)=beta + if (standardMap->getLocalElement(standardMap->getMaxGlobalIndex()) >= 0) { + int lastDof = standardMap->getLocalElement(standardMap->getMaxGlobalIndex()); + if ((flag == F_ONLY) || (flag == ALL)) + (*rhs)[lastDof] = (*soln)[lastDof] - beta; + if ((flag == MATRIX_ONLY) || (flag == ALL)) { + int row = standardMap->getMaxGlobalIndex(); + column = row; + jac = 1.0 * jac_coeff; + A->replaceGlobalValues(row, 1, &jac, &column); + column = row - 1; + jac = 0.0 * jac_coeff; + A->replaceGlobalValues(row, 1, &jac, &column); + } + } + // Sync up processors to be safe + comm->barrier(); + + A->fillComplete(); + + return true; +} + +Tpetra_Vector& Pitchfork_FiniteElementProblem::getSolution() { return *initialSolution; } + +Tpetra_CrsMatrix& Pitchfork_FiniteElementProblem::getJacobian() { return *A; } + +bool Pitchfork_FiniteElementProblem::setParameter(std::string label, double value) { + if (label == "lambda") + lambda = value; + else if (label == "alpha") + alpha = value; + else if (label == "beta") + beta = value; + else if (label == "Homotopy Continuation Parameter") { + // do nothing for now + } else { + // do nothing (may be a constraint parameter that we don't know about) + } + return true; +} + +tcsrgraph_t& Pitchfork_FiniteElementProblem::generateGraph(tcsrgraph_t& AAA) { + // Declare required variables + int i, j; + int row, column; + int overlapNumLocalElements = overlapMap->NumMyElements(); + + // Loop Over # of Finite Elements on Processor + for (int ne = 0; ne < overlapNumLocalElements - 1; ne++) { + // Loop over Nodes in Element + for (i = 0; i < 2; i++) { + row = overlapMap->GID(ne + i); + + // Loop over Trial Functions + for (j = 0; j < 2; j++) { + // If this row is owned by current processor, add the index + if (standardMap->MyGID(row)) { + column = overlapMap->GID(ne + j); + AAA.insertGlobalIndices(row, 1, &column); + } + } + } + } + AAA.FillComplete(); + return AAA; +} + +double Pitchfork_FiniteElementProblem::source_term(double x) { + return lambda * x - alpha * x * x + beta * x * x * x; +} + +double Pitchfork_FiniteElementProblem::source_deriv(double x) { + return lambda - 2.0 * alpha * x + 3.0 * beta * x * x; +} diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/Pitchfork_FiniteElementProblem.hpp b/packages/nox/test/tpetra/LOCA_TestProblems/Pitchfork_FiniteElementProblem.hpp new file mode 100644 index 000000000000..94d7a4564c0a --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/Pitchfork_FiniteElementProblem.hpp @@ -0,0 +1,153 @@ +// $Id$ +// $Source$ + +//@HEADER +// ************************************************************************ +// +// LOCA: Library of Continuation Algorithms Package +// Copyright (2005) 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 + +// Finite Element Problem Class +/* Provides function (F) and Jacobian evaluations for the following equation + * via a 1D linear finite element discretization with Epetra objects. + * + * d2u + * --- + lambda * u - alpha * u**2 + beta * u**3 = 0 + * dx2 + * + * subject to the boundar conditions u(-1) = u(1) = beta. + */ + +#ifndef _PITCHFORK_FINITEELEMENTPROBLEM_H +#define _PITCHFORK_FINITEELEMENTPROBLEM_H + +#include "FiniteElementProblem.hpp" // base class + +// Forward Declarations +namespace Teuchos { template class Comm; } + +// Finite Element Problem Class +template +class Pitchfork_FiniteElementProblem : public FiniteElementProblem { + +public: + + // Types + using scalar_type = typename FiniteElementProblem::scalar_type; + using vector_type = typename FiniteElementProblem::vector_type; + using row_matrix_type = typename FiniteElementProblem::row_matrix_type; + using crs_matrix_type = typename FiniteElementProblem::crs_matrix_type; + using import_type = typename FiniteElementProblem::import_type; + using csr_graph_type = typename FiniteElementProblem::csr_graph_type; + using map_type = typename FiniteElementProblem::map_type; + + // Constructor + Pitchfork_FiniteElementProblem(int numGlobalElements, Teuchos::RCP>& comm); + + // Destructor + virtual ~Pitchfork_FiniteElementProblem(); + + // Evaluates the function (F) and/or the Jacobian using the solution + // values in solnVector. + virtual bool evaluate(FillType f, const vector_type *solnVector, + vector_type *rhsVector, row_matrix_type *matrix, + scalar_type jac_coeff = 1.0, + scalar_type mass_coeff = 0.0); + + // Return a reference to the Tpetra_Vector with the initial guess + // that is generated by the FiniteElementProblem class. + virtual vector_type& getSolution(); + + // Return a reference to the Tpetra_Vector with the Jacobian + // that is generated by the FiniteElementProblem class. + virtual crs_matrix_type & getJacobian(); + + // Set a bifurcation parameter in the application physics + virtual bool setParameter(std::string label, scalar_type value); + +private: + + // inserts the global column indices into the Graph + csr_graph_type& generateGraph(csr_graph_type& AA); + + // Computes the source term + scalar_type source_term(scalar_type x); + + // Computes the derivative of the source term + scalar_type source_deriv(scalar_type x); + + // Private to prohibit copying + Pitchfork_FiniteElementProblem(const Pitchfork_FiniteElementProblem&); + + // Private to prohibit copying + Pitchfork_FiniteElementProblem& operator=(const Pitchfork_FiniteElementProblem&); + +private: + + FillType flag; + map_type *standardMap; + map_type *overlapMap; + import_type *importer; + vector_type *initialSolution; + vector_type *rhs; + csr_graph_type *AA; + crs_matrix_type *A; + Teuchos::RCP> comm; + + int rank; // Process number + int numProcs; // Total number of processes + int numLocalElements; // Number of elements owned by this process + int numGlobalElements; // Total Number of elements + + public: + + scalar_type lambda; // factor used on nonlinear term + scalar_type alpha; // factor used on nonlinear term + scalar_type beta; // factor used on nonlinear term +}; +#endif + + + + diff --git a/packages/nox/test/tpetra/LOCA_UnitTests/CMakeLists.txt b/packages/nox/test/tpetra/LOCA_UnitTests/CMakeLists.txt new file mode 100644 index 000000000000..64ff032ccf36 --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_UnitTests/CMakeLists.txt @@ -0,0 +1,90 @@ + + +TRIBITS_INCLUDE_DIRECTORIES(REQUIRED_DURING_INSTALLATION_TESTING ${CMAKE_CURRENT_SOURCE_DIR}) + +TRIBITS_INCLUDE_DIRECTORIES(REQUIRED_DURING_INSTALLATION_TESTING ${CMAKE_CURRENT_SOURCE_DIR}/../LOCA_TestProblems) + +TRIBITS_INCLUDE_DIRECTORIES(REQUIRED_DURING_INSTALLATION_TESTING ${CMAKE_CURRENT_SOURCE_DIR}/../../utils) + +IF(NOX_ENABLE_LOCA) # IF(NOX_ENABLE_ABSTRACT_IMPLEMENTATION_EPETRA AND NOX_ENABLE_LOCA) + + # NGA Dev Tests (to be removed as soon as locatpetratestproblems library is completed) + TRIBITS_ADD_EXECUTABLE_AND_TEST( + Tpetra_LOCA_NGA_DevTests + SOURCES + NGA_DevTests.cpp + COMM serial mpi + NUM_MPI_PROCS 2 + TESTONLYLIBS noxtestutils locatpetratestproblems + ARGS -v + STANDARD_PASS_OUTPUT + ) + +# TRIBITS_ADD_EXECUTABLE_AND_TEST( +# LOCA_HouseholderBorderedSolve +# SOURCES +# HouseholderBorderedSolve.C +# COMM serial mpi +# NUM_MPI_PROCS 2 +# TESTONLYLIBS noxtestutils locaepetratestproblems +# ARGS -v +# PASS_REGULAR_EXPRESSION "All tests passed" +# ) + +# TRIBITS_ADD_EXECUTABLE_AND_TEST( +# LOCA_HouseholderTransposeBorderedSolve +# SOURCES +# HouseholderTransposeBorderedSolve.C +# COMM serial mpi +# NUM_MPI_PROCS 2 +# TESTONLYLIBS noxtestutils locaepetratestproblems +# ARGS -v +# PASS_REGULAR_EXPRESSION "All tests passed" +# ) + +# TRIBITS_ADD_EXECUTABLE_AND_TEST( +# LOCA_TransposeSolve +# SOURCES +# TransposeSolve.C +# COMM serial mpi +# NUM_MPI_PROCS 2 +# TESTONLYLIBS noxtestutils locaepetratestproblems +# ARGS -v +# PASS_REGULAR_EXPRESSION "All tests passed" +# ) + + # WIP + # TRIBITS_ADD_EXECUTABLE_AND_TEST( + # Tpetra_LOCA_CompositeConstraint + # SOURCES + # CompositeConstraint.C + # COMM serial mpi + # NUM_MPI_PROCS 2 + # TESTONLYLIBS noxtestutils locatpetratestproblems + # ARGS -v + # STANDARD_PASS_OUTPUT + # ) + +# TRIBITS_ADD_EXECUTABLE_AND_TEST( +# LOCA_CompositeConstraintMVDX +# SOURCES +# CompositeConstraintMVDX.C +# COMM serial mpi +# NUM_MPI_PROCS 2 +# TESTONLYLIBS noxtestutils locaepetratestproblems +# ARGS -v +# PASS_REGULAR_EXPRESSION "All tests passed" +# ) + +# TRIBITS_ADD_EXECUTABLE_AND_TEST( +# LOCA_NaturalContResidualFills +# SOURCES +# NaturalContResidualFills.C +# COMM serial mpi +# NUM_MPI_PROCS 2 +# TESTONLYLIBS noxtestutils locaepetratestproblems +# ARGS -v +# PASS_REGULAR_EXPRESSION "All tests passed" +# ) + +ENDIF() diff --git a/packages/nox/test/tpetra/LOCA_UnitTests/NGA_DevTests.cpp b/packages/nox/test/tpetra/LOCA_UnitTests/NGA_DevTests.cpp new file mode 100644 index 000000000000..95fbbb5fd7a6 --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_UnitTests/NGA_DevTests.cpp @@ -0,0 +1,29 @@ +// Teuchos +#include +#include + +// Tpetra +#include + +// Lib headers +#include "Pitchfork_FiniteElementProblem.hpp" +#include "LOCALinearConstraint.hpp" + +// This is a temporary test created to run any code from the test libs +// to checkk if it is running +int main(int argc, char* argv[]) { + + using ST = typename Tpetra::MultiVector::scalar_type; + + Teuchos::GlobalMPISession mpiSession(&argc, &argv); + auto comm = Tpetra::getDefaultComm(); + const int myRank = comm->getRank(); + + // test Pitchfork_FiniteElementProblem (inheriting FiniteElementProblem) + int numGlobalElements = 1000; + auto feProblem = Pitchfork_FiniteElementProblem(numGlobalElements, comm); + + if (myRank == 0) { + std::cout << "End Result: TEST PASSED" << std::endl; + } +} \ No newline at end of file