From 062f0b5916d6e1f2aa513988609309a4251a341a Mon Sep 17 00:00:00 2001 From: antoinemeyer5 Date: Mon, 18 Sep 2023 09:21:23 -0400 Subject: [PATCH 01/13] #197: Init addition --- packages/nox/test/tpetra/CMakeLists.txt | 2 + .../test/tpetra/LOCA_TestProblems/Basis.cpp | 101 ++++++++++++++++++ .../test/tpetra/LOCA_TestProblems/Basis.hpp | 81 ++++++++++++++ .../tpetra/LOCA_TestProblems/CMakeLists.txt | 36 +++++++ 4 files changed, 220 insertions(+) create mode 100644 packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp create mode 100644 packages/nox/test/tpetra/LOCA_TestProblems/Basis.hpp create mode 100644 packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt diff --git a/packages/nox/test/tpetra/CMakeLists.txt b/packages/nox/test/tpetra/CMakeLists.txt index 5334a5313c19..a79074ec3c02 100644 --- a/packages/nox/test/tpetra/CMakeLists.txt +++ b/packages/nox/test/tpetra/CMakeLists.txt @@ -68,3 +68,5 @@ IF(NOX_ENABLE_ABSTRACT_IMPLEMENTATION_THYRA AND ENDIF() ENDIF() + +ADD_SUBDIRECTORY(LOCA_TestProblems) 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..a532fb55277a --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp @@ -0,0 +1,101 @@ +// $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 "NOX_Common.H" +#include "Basis.H" +#include +#include + +// Constructor +Basis::Basis() : + phi(NULL), + dphide(NULL), + uu(0.0), + xx(0.0), + duu(0.0), + eta(0.0), + wt(0.0), + dx(0.0) +{ + phi = new double[2]; + dphide = new double[2]; +} + +// Destructor +Basis::~Basis() { + delete [] phi; + delete [] dphide; +} + +// Calculates a linear 1D basis +void Basis::getBasis(int gp, double *x, double *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..2a2c31a986b4 --- /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 + +class Basis { + + public: + + // Constructor + Basis(); + + // Destructor + ~Basis(); + + // Calculates the values of u and x at the specified gauss point + void getBasis(int gp, double *x, double *u); + + private: + // Private to prohibit copying + Basis(const Basis&); + Basis& operator=(const Basis&); + + public: + // Variables that are calculated at the gauss point + double *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..52df50000ec2 --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt @@ -0,0 +1,36 @@ + +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 + NormConstraint.hpp + LOCALinearConstraint.hpp +) + +APPEND_SET(SOURCES + Basis.cpp + Problem_Interface.cpp + Tcubed_FiniteElementProblem.cpp + Pitchfork_FiniteElementProblem.cpp + NormConstraint.cpp + LOCALinearConstraint.cpp +) + +IF(NOX_ENABLE_ABSTRACT_IMPLEMENTATION_TPETRA AND NOX_ENABLE_LOCA) + + TRIBITS_ADD_LIBRARY( + locatpetratestproblems + HEADERS ${HEADERS} + SOURCES ${SOURCES} + TESTONLY + DEPLIBS loca locatpetra + ) + +ENDIF() From 3e83affb75a06e89873d1836afa31df6ba4f3793 Mon Sep 17 00:00:00 2001 From: antoinemeyer5 Date: Mon, 18 Sep 2023 10:18:29 -0400 Subject: [PATCH 02/13] #197: Build ok with basis class --- .../test/tpetra/LOCA_TestProblems/Basis.cpp | 4 ++-- .../tpetra/LOCA_TestProblems/CMakeLists.txt | 24 +++++++++---------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp b/packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp index a532fb55277a..5801b0007473 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp +++ b/packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp @@ -48,8 +48,8 @@ // ************************************************************************ //@HEADER -#include "NOX_Common.H" -#include "Basis.H" +// #include "NOX_Common.H" +#include "Basis.hpp" #include #include diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt index 52df50000ec2..84842a082693 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt +++ b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt @@ -6,24 +6,24 @@ TRIBITS_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}) APPEND_SET(HEADERS Basis.hpp - FiniteElementProblem.hpp - Problem_Interface.hpp - Tcubed_FiniteElementProblem.hpp - Pitchfork_FiniteElementProblem.hpp - NormConstraint.hpp - LOCALinearConstraint.hpp +# FiniteElementProblem.hpp +# Problem_Interface.hpp +# Tcubed_FiniteElementProblem.hpp +# Pitchfork_FiniteElementProblem.hpp +# NormConstraint.hpp +# LOCALinearConstraint.hpp ) APPEND_SET(SOURCES Basis.cpp - Problem_Interface.cpp - Tcubed_FiniteElementProblem.cpp - Pitchfork_FiniteElementProblem.cpp - NormConstraint.cpp - LOCALinearConstraint.cpp +# Problem_Interface.cpp +# Tcubed_FiniteElementProblem.cpp +# Pitchfork_FiniteElementProblem.cpp +# NormConstraint.cpp +# LOCALinearConstraint.cpp ) -IF(NOX_ENABLE_ABSTRACT_IMPLEMENTATION_TPETRA AND NOX_ENABLE_LOCA) +IF(#[[Change this `NOX_ENABLE_ABSTRACT_IMPLEMENTATION_EPETRA` for Tpetra]] #[[AND]] NOX_ENABLE_LOCA) TRIBITS_ADD_LIBRARY( locatpetratestproblems From f6fcfef55682a4eaa853ff3ccff854883ef617f4 Mon Sep 17 00:00:00 2001 From: antoinemeyer5 Date: Mon, 18 Sep 2023 10:55:41 -0400 Subject: [PATCH 03/13] #197: Template the Basis class --- .../test/tpetra/LOCA_TestProblems/Basis.cpp | 101 ------------------ .../test/tpetra/LOCA_TestProblems/Basis.hpp | 74 ++++++++++--- 2 files changed, 57 insertions(+), 118 deletions(-) delete mode 100644 packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp b/packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp deleted file mode 100644 index 5801b0007473..000000000000 --- a/packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp +++ /dev/null @@ -1,101 +0,0 @@ -// $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 "NOX_Common.H" -#include "Basis.hpp" -#include -#include - -// Constructor -Basis::Basis() : - phi(NULL), - dphide(NULL), - uu(0.0), - xx(0.0), - duu(0.0), - eta(0.0), - wt(0.0), - dx(0.0) -{ - phi = new double[2]; - dphide = new double[2]; -} - -// Destructor -Basis::~Basis() { - delete [] phi; - delete [] dphide; -} - -// Calculates a linear 1D basis -void Basis::getBasis(int gp, double *x, double *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 index 2a2c31a986b4..ad15f080dd67 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/Basis.hpp +++ b/packages/nox/test/tpetra/LOCA_TestProblems/Basis.hpp @@ -53,29 +53,69 @@ #ifndef _NOX_EXAMPLE_TPETRA_LINEAR_BASIS_H #define _NOX_EXAMPLE_TPETRA_LINEAR_BASIS_H -class Basis { +// #include "NOX_Common.H" +#include +#include +#include - public: +template +class Basis +{ + using stdVector = typename std::vector; - // Constructor - Basis(); + public: + // Constructor + Basis() : + phi(NULL), + dphide(NULL), + uu(0.0), + xx(0.0), + duu(0.0), + eta(0.0), + wt(0.0), + dx(0.0) + { + phi(2); + dphide(2); + } - // Destructor - ~Basis(); + // Calculates the values of u and x at the specified gauss point + void getBasis(int gp, stdVector x, stdVector 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;} - // Calculates the values of u and x at the specified gauss point - void getBasis(int gp, double *x, double *u); + // 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; - private: - // Private to prohibit copying - Basis(const Basis&); - Basis& operator=(const Basis&); + // 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]; + } - public: - // Variables that are calculated at the gauss point - double *phi, *dphide; - double uu, xx, duu, eta, wt; - double dx; + return; + } + + private: + // Private to prohibit copying + Basis(const Basis&); + Basis& operator=(const Basis&); + + public: + // Variables that are calculated at the gauss point + stdVector phi, dphide; + ScalarType uu, xx, duu, eta, wt; + ScalarType dx; }; #endif From e0b8a6865922824476232cf55f7d927c7d13714e Mon Sep 17 00:00:00 2001 From: antoinemeyer5 Date: Mon, 18 Sep 2023 10:55:56 -0400 Subject: [PATCH 04/13] #197: Template the Basis class --- packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt index 84842a082693..8ca8427b8c0b 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt +++ b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt @@ -15,7 +15,6 @@ APPEND_SET(HEADERS ) APPEND_SET(SOURCES - Basis.cpp # Problem_Interface.cpp # Tcubed_FiniteElementProblem.cpp # Pitchfork_FiniteElementProblem.cpp From c16cd1401992eacfe6c8951d57c34daafce8e061 Mon Sep 17 00:00:00 2001 From: antoinemeyer5 Date: Mon, 18 Sep 2023 11:36:12 -0400 Subject: [PATCH 05/13] #197: Undo template the Basis class because of the library --- .../test/tpetra/LOCA_TestProblems/Basis.cpp | 89 +++++++++++++++++++ .../test/tpetra/LOCA_TestProblems/Basis.hpp | 52 ++--------- .../tpetra/LOCA_TestProblems/CMakeLists.txt | 1 + 3 files changed, 96 insertions(+), 46 deletions(-) create mode 100644 packages/nox/test/tpetra/LOCA_TestProblems/Basis.cpp 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 index ad15f080dd67..1f76a0a8f5cf 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/Basis.hpp +++ b/packages/nox/test/tpetra/LOCA_TestProblems/Basis.hpp @@ -53,58 +53,18 @@ #ifndef _NOX_EXAMPLE_TPETRA_LINEAR_BASIS_H #define _NOX_EXAMPLE_TPETRA_LINEAR_BASIS_H -// #include "NOX_Common.H" +#include #include #include -#include -template class Basis { - using stdVector = typename std::vector; - public: // Constructor - Basis() : - phi(NULL), - dphide(NULL), - uu(0.0), - xx(0.0), - duu(0.0), - eta(0.0), - wt(0.0), - dx(0.0) - { - phi(2); - dphide(2); - } + Basis(); // Calculates the values of u and x at the specified gauss point - void getBasis(int gp, stdVector x, stdVector 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; - } + void getBasis(int gp, std::vector x, std::vector u); private: // Private to prohibit copying @@ -113,9 +73,9 @@ class Basis public: // Variables that are calculated at the gauss point - stdVector phi, dphide; - ScalarType uu, xx, duu, eta, wt; - ScalarType dx; + 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 index 8ca8427b8c0b..84842a082693 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt +++ b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt @@ -15,6 +15,7 @@ APPEND_SET(HEADERS ) APPEND_SET(SOURCES + Basis.cpp # Problem_Interface.cpp # Tcubed_FiniteElementProblem.cpp # Pitchfork_FiniteElementProblem.cpp From 4b228e1bc1ffa6570e08e969228ae26450c58099 Mon Sep 17 00:00:00 2001 From: antoinemeyer5 Date: Tue, 19 Sep 2023 19:36:25 +0200 Subject: [PATCH 06/13] #197: Add FiniteElementProblem class --- .../FiniteElementProblem.cpp | 337 ++++++++++++++++++ 1 file changed, 337 insertions(+) create mode 100644 packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp 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..b42777470586 --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp @@ -0,0 +1,337 @@ +//@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 "NOX_Common.H" +#include "Epetra_Comm.h" +#include "Epetra_Map.h" +#include "Epetra_Vector.h" +#include "Epetra_Import.h" +#include "Epetra_CrsGraph.h" +#include "Epetra_CrsMatrix.h" +#include "Basis.H" + +#include "FiniteElementProblem.H" + +// Constructor - creates the Epetra objects (maps and vectors) +FiniteElementProblem::FiniteElementProblem(int numGlobalElements, + Epetra_Comm& comm, + double s) : + Comm(&comm), + NumGlobalElements(numGlobalElements), + scale(s) +{ + + // Commonly used variables + int i; + MyPID = Comm->MyPID(); // Process ID + NumProc = 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 Epetra_Map(NumGlobalElements, 0, *Comm); + + // Get the number of elements owned by this processor + NumMyElements = StandardMap->NumMyElements(); + + // 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 Epetra_Map(*StandardMap); + } else { + + int OverlapNumMyElements; + int OverlapMinMyGID; + OverlapNumMyElements = NumMyElements + 2; + if ((MyPID == 0) || (MyPID == NumProc - 1)) + OverlapNumMyElements --; + + if (MyPID==0) + OverlapMinMyGID = StandardMap->MinMyGID(); + else + OverlapMinMyGID = StandardMap->MinMyGID() - 1; + + int* OverlapMyGlobalElements = new int[OverlapNumMyElements]; + + for (i = 0; i < OverlapNumMyElements; i ++) + OverlapMyGlobalElements[i] = OverlapMinMyGID + i; + + OverlapMap = new Epetra_Map(-1, OverlapNumMyElements, + OverlapMyGlobalElements, 0, *Comm); + + delete [] OverlapMyGlobalElements; + + } // End Overlap map construction ************************************* + + // Construct Linear Objects + Importer = new Epetra_Import(*OverlapMap, *StandardMap); + initialSolution = new Epetra_Vector(*StandardMap); + AA = new Epetra_CrsGraph(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 Epetra_CrsMatrix (Copy, *AA); + A->FillComplete(); + + // Set default bifurcation values + factor = 0.1; + leftBC = 1.0; + rightBC = 1.0; +} + +// Destructor +FiniteElementProblem::~FiniteElementProblem() +{ + delete AA; + delete A; + delete initialSolution; + delete Importer; + delete OverlapMap; + delete StandardMap; +} + +// Matrix and Residual Fills +bool FiniteElementProblem::evaluate(FillType f, + const Epetra_Vector* soln, + Epetra_Vector* tmp_rhs, + Epetra_RowMatrix* 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 + Epetra_Vector u(*OverlapMap); + Epetra_Vector x(*OverlapMap); + + // Export Solution to Overlap vector + u.Import(*soln, *Importer, Insert); + + // Declare required variables + int i,j,ierr; + int OverlapNumMyElements = OverlapMap->NumMyElements(); + + int OverlapMinMyGID; + if (MyPID==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=1.0; + double dx=Length/((double) NumGlobalElements-1); + for (i=0; i < OverlapNumMyElements; i++) { + x[i]=dx*((double) 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->GID(ne+i); + //printf("Proc=%d GlobalRow=%d LocalRow=%d Owned=%d\n", + // MyPID, 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->MyGID(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 (MyPID==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->MaxAllGID()) >= 0 ) { + int lastDof = StandardMap->LID(StandardMap->MaxAllGID()); + if ((flag == F_ONLY) || (flag == ALL)) + (*rhs)[lastDof] = (*soln)[lastDof] - rightBC; + if ((flag == MATRIX_ONLY) || (flag == ALL)) { + int row = StandardMap->MaxAllGID(); + 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; +} + +Epetra_Vector& FiniteElementProblem::getSolution() +{ + return *initialSolution; +} + +Epetra_CrsMatrix& FiniteElementProblem::getJacobian() +{ + return *A; +} + +bool FiniteElementProblem::setParameter(std::string label, double 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; +} + +Epetra_CrsGraph& FiniteElementProblem::generateGraph(Epetra_CrsGraph& AAA) +{ + + // Declare required variables + int i,j; + int row, column; + int OverlapNumMyElements = OverlapMap->NumMyElements(); + int OverlapMinMyGID; + if (MyPID==0) OverlapMinMyGID = StandardMap->MinMyGID(); + else OverlapMinMyGID = StandardMap->MinMyGID()-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(); +// AAA.SortIndices(); +// AAA.RemoveRedundantIndices(); + return AAA; +} +*/ From 5abaf2d6462ff12b2b43438ad9d5e0533659293a Mon Sep 17 00:00:00 2001 From: antoinemeyer5 Date: Tue, 19 Sep 2023 19:45:30 +0200 Subject: [PATCH 07/13] #197: Add FiniteElementProblem class - wip --- .../nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp index b42777470586..9c684206d2c3 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp +++ b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp @@ -334,4 +334,5 @@ Epetra_CrsGraph& FiniteElementProblem::generateGraph(Epetra_CrsGraph& AAA) // AAA.RemoveRedundantIndices(); return AAA; } + */ From aa20662424e058875b3f75483cece2823dc95db8 Mon Sep 17 00:00:00 2001 From: antoinemeyer5 Date: Wed, 20 Sep 2023 11:21:17 +0200 Subject: [PATCH 08/13] #197: Add FiniteElementProblem class with template --- .../tpetra/LOCA_TestProblems/CMakeLists.txt | 7 +- .../FiniteElementProblem.cpp | 157 ++++++++---------- .../FiniteElementProblem.hpp | 106 ++++++++++++ 3 files changed, 180 insertions(+), 90 deletions(-) create mode 100644 packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.hpp diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt index 84842a082693..6f2cb03f4396 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt +++ b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt @@ -6,7 +6,7 @@ TRIBITS_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}) APPEND_SET(HEADERS Basis.hpp -# FiniteElementProblem.hpp + FiniteElementProblem.hpp # Problem_Interface.hpp # Tcubed_FiniteElementProblem.hpp # Pitchfork_FiniteElementProblem.hpp @@ -23,7 +23,10 @@ APPEND_SET(SOURCES # LOCALinearConstraint.cpp ) -IF(#[[Change this `NOX_ENABLE_ABSTRACT_IMPLEMENTATION_EPETRA` for Tpetra]] #[[AND]] NOX_ENABLE_LOCA) +# Add this below: +# IF(NOX_ENABLE_ABSTRACT_IMPLEMENTATION_EPETRA AND NOX_ENABLE_LOCA) + +IF(NOX_ENABLE_LOCA) TRIBITS_ADD_LIBRARY( locatpetratestproblems diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp index 9c684206d2c3..4073cc51ba22 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp +++ b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp @@ -45,22 +45,16 @@ // ************************************************************************ //@HEADER -/* -#include "NOX_Common.H" -#include "Epetra_Comm.h" -#include "Epetra_Map.h" -#include "Epetra_Vector.h" -#include "Epetra_Import.h" -#include "Epetra_CrsGraph.h" -#include "Epetra_CrsMatrix.h" -#include "Basis.H" - -#include "FiniteElementProblem.H" - -// Constructor - creates the Epetra objects (maps and vectors) -FiniteElementProblem::FiniteElementProblem(int numGlobalElements, - Epetra_Comm& comm, - double s) : +#include + +#include "Basis.hpp" +#include "FiniteElementProblem.hpp" + +// Constructor - creates the Tpetra objects (maps and vectors) +template typename +FiniteElementProblem::FiniteElementProblem(int numGlobalElements, + Teuchos::RCP> &comm, + ST s) : Comm(&comm), NumGlobalElements(numGlobalElements), scale(s) @@ -68,20 +62,20 @@ FiniteElementProblem::FiniteElementProblem(int numGlobalElements, // Commonly used variables int i; - MyPID = Comm->MyPID(); // Process ID - NumProc = Comm->NumProc(); // Total number of processes + MyPID = 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 Epetra_Map(NumGlobalElements, 0, *Comm); + StandardMap = new tmap_t(NumGlobalElements, 0, *Comm); // Get the number of elements owned by this processor - NumMyElements = StandardMap->NumMyElements(); + 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 Epetra_Map(*StandardMap); + OverlapMap = new tmap_t(*StandardMap); } else { int OverlapNumMyElements; @@ -91,26 +85,24 @@ FiniteElementProblem::FiniteElementProblem(int numGlobalElements, OverlapNumMyElements --; if (MyPID==0) - OverlapMinMyGID = StandardMap->MinMyGID(); + OverlapMinMyGID = StandardMap->getMinGlobalIndex(); else - OverlapMinMyGID = StandardMap->MinMyGID() - 1; + OverlapMinMyGID = StandardMap->getMinGlobalIndex() - 1; int* OverlapMyGlobalElements = new int[OverlapNumMyElements]; for (i = 0; i < OverlapNumMyElements; i ++) OverlapMyGlobalElements[i] = OverlapMinMyGID + i; - OverlapMap = new Epetra_Map(-1, OverlapNumMyElements, + OverlapMap = new tmap_t(-1, OverlapNumMyElements, OverlapMyGlobalElements, 0, *Comm); - delete [] OverlapMyGlobalElements; - } // End Overlap map construction ************************************* // Construct Linear Objects - Importer = new Epetra_Import(*OverlapMap, *StandardMap); - initialSolution = new Epetra_Vector(*StandardMap); - AA = new Epetra_CrsGraph(Copy, *StandardMap, 5); + 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); @@ -118,7 +110,7 @@ FiniteElementProblem::FiniteElementProblem(int numGlobalElements, // 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 Epetra_CrsMatrix (Copy, *AA); + A = new tcrsmatrix_t(Copy, *AA); A->FillComplete(); // Set default bifurcation values @@ -127,22 +119,12 @@ FiniteElementProblem::FiniteElementProblem(int numGlobalElements, rightBC = 1.0; } -// Destructor -FiniteElementProblem::~FiniteElementProblem() -{ - delete AA; - delete A; - delete initialSolution; - delete Importer; - delete OverlapMap; - delete StandardMap; -} - // Matrix and Residual Fills +template typename bool FiniteElementProblem::evaluate(FillType f, - const Epetra_Vector* soln, - Epetra_Vector* tmp_rhs, - Epetra_RowMatrix* tmp_matrix) + const tvector_t* soln, + tvector_t* tmp_rhs, + trowmatrix_t* tmp_matrix) { flag = f; @@ -150,41 +132,41 @@ bool FiniteElementProblem::evaluate(FillType f, if (flag == F_ONLY) { rhs = tmp_rhs; } else if (flag == MATRIX_ONLY) { - A = dynamic_cast (tmp_matrix); + A = dynamic_cast (tmp_matrix); } else if (flag == ALL) { rhs = tmp_rhs; - A = dynamic_cast (tmp_matrix); + 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 - Epetra_Vector u(*OverlapMap); - Epetra_Vector x(*OverlapMap); + tvector_t u(*OverlapMap); + tvector_t x(*OverlapMap); // Export Solution to Overlap vector - u.Import(*soln, *Importer, Insert); + u.doImport(*soln, *Importer, Tpetra::INSERT); // Declare required variables int i,j,ierr; int OverlapNumMyElements = OverlapMap->NumMyElements(); int OverlapMinMyGID; - if (MyPID==0) OverlapMinMyGID = StandardMap->MinMyGID(); - else OverlapMinMyGID = StandardMap->MinMyGID()-1; + if (MyPID==0) OverlapMinMyGID = StandardMap->getMinGlobalIndex(); + else OverlapMinMyGID = StandardMap->getMinGlobalIndex()-1; int row, column; - double jac; - double xx[2]; - double uu[2]; + ST jac; + ST xx[2]; + ST uu[2]; Basis basis; // Create the nodal coordinates - double Length=1.0; - double dx=Length/((double) NumGlobalElements-1); + ST Length=1.0; + ST dx=Length/((ST) NumGlobalElements-1); for (i=0; i < OverlapNumMyElements; i++) { - x[i]=dx*((double) OverlapMinMyGID+i); + x[i]=dx*((ST) OverlapMinMyGID+i); } // Zero out the objects that will be filled @@ -226,7 +208,7 @@ bool FiniteElementProblem::evaluate(FillType f, 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); + ierr=A->sumIntoGlobalValues(row, 1, &jac, &column); } } } @@ -242,46 +224,49 @@ bool FiniteElementProblem::evaluate(FillType f, if ((flag == MATRIX_ONLY) || (flag == ALL)) { column=0; jac=1.0; - A->ReplaceGlobalValues(0, 1, &jac, &column); + A->replaceGlobalValues(0, 1, &jac, &column); column=1; jac=0.0; - A->ReplaceGlobalValues(0, 1, &jac, &column); + A->replaceGlobalValues(0, 1, &jac, &column); } } - if ( StandardMap->LID(StandardMap->MaxAllGID()) >= 0 ) { - int lastDof = StandardMap->LID(StandardMap->MaxAllGID()); + if ( StandardMap->LID(StandardMap->getMaxAllGlobalIndex()) >= 0 ) { + int lastDof = StandardMap->LID(StandardMap->getMaxAllGlobalIndex()); if ((flag == F_ONLY) || (flag == ALL)) (*rhs)[lastDof] = (*soln)[lastDof] - rightBC; if ((flag == MATRIX_ONLY) || (flag == ALL)) { - int row = StandardMap->MaxAllGID(); + int row = StandardMap->getMaxAllGlobalIndex(); column = row; jac=1.0; - A->ReplaceGlobalValues(row, 1, &jac, &column); + A->replaceGlobalValues(row, 1, &jac, &column); column=row-1; jac=0.0; - A->ReplaceGlobalValues(row, 1, &jac, &column); + A->replaceGlobalValues(row, 1, &jac, &column); } } // Sync up processors to be safe - Comm->Barrier(); + Comm->barrier(); - A->FillComplete(); + A->fillComplete(); return true; } -Epetra_Vector& FiniteElementProblem::getSolution() +template typename +Teuchos::RCP FiniteElementProblem::getSolution() { - return *initialSolution; + return initialSolution; } -Epetra_CrsMatrix& FiniteElementProblem::getJacobian() +template typename +Teuchos::RCP FiniteElementProblem::getJacobian() { - return *A; + return A; } -bool FiniteElementProblem::setParameter(std::string label, double value) +template typename +bool FiniteElementProblem::setParameter(std::string label, ST value) { if (label == "Nonlinear Factor") factor = value; @@ -300,16 +285,16 @@ bool FiniteElementProblem::setParameter(std::string label, double value) return true; } -Epetra_CrsGraph& FiniteElementProblem::generateGraph(Epetra_CrsGraph& AAA) +template typename +tcrsgraph_t& FiniteElementProblem::generateGraph(tcrsgraph_t& AAA) { - // Declare required variables int i,j; int row, column; - int OverlapNumMyElements = OverlapMap->NumMyElements(); + int OverlapNumMyElements = OverlapMap->getLocalNumElements(); int OverlapMinMyGID; - if (MyPID==0) OverlapMinMyGID = StandardMap->MinMyGID(); - else OverlapMinMyGID = StandardMap->MinMyGID()-1; + if (MyPID==0) OverlapMinMyGID = StandardMap->getMinGlobalIndex(); + else OverlapMinMyGID = StandardMap->getMinGlobalIndex()-1; // Loop Over # of Finite Elements on Processor for (int ne=0; ne < OverlapNumMyElements-1; ne++) { @@ -320,19 +305,15 @@ Epetra_CrsGraph& FiniteElementProblem::generateGraph(Epetra_CrsGraph& AAA) // 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); - } + + // 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(); -// AAA.SortIndices(); -// AAA.RemoveRedundantIndices(); + 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..70c0378cd437 --- /dev/null +++ b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.hpp @@ -0,0 +1,106 @@ +// $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 ST = typename Tpetra::Vector::scalar_type; + using LO = typename Tpetra::Vector<>::local_ordinal_type; + using GO = typename Tpetra::Vector<>::global_ordinal_type; + using NO = typename Tpetra::Map<>::node_type; + + using tvector_t = typename Tpetra::Vector; + using trowmatrix_t = typename Tpetra::RowMatrix; + using tcrsmatrix_t = typename Tpetra::CrsMatrix; + +public: + + // Constructor + FiniteElementProblem() {} + + // Evaluates the function (F) and/or the Jacobian using the solution + // values in solnVector. + virtual bool evaluate(FillType f, const tvector_t *solnVector, + tvector_t *rhsVector, + trowmatrix_t *matrix, + ST jac_coeff = 1.0, + ST mass_coeff = 0.0) = 0; + + // Return a reference to the Tpetra_Vector with the initial guess + // that is generated by the FiniteElementProblem class. + virtual tvector_t& getSolution() = 0; + + // Return a reference to the Tpetra_Vector with the Jacobian + // that is generated by the FiniteElementProblem class. + virtual tcrsmatrix_t& getJacobian() = 0; + + // Set a bifurcation parameter in the application physics + virtual bool setParameter(std::string label, ST value) = 0; + +}; +#endif From b51539f7af48cfe5be4e01c1637fcdfb79451ded Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20Dutheillet-Lamonth=C3=A9zie?= Date: Fri, 29 Sep 2023 16:25:40 +0200 Subject: [PATCH 09/13] #197: update FiniteElementProblem class --- .../FiniteElementProblem.cpp | 257 +++++++++--------- .../FiniteElementProblem.hpp | 31 ++- 2 files changed, 141 insertions(+), 147 deletions(-) diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp index 4073cc51ba22..a2ac87da976b 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp +++ b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.cpp @@ -48,61 +48,58 @@ #include #include "Basis.hpp" + #include "FiniteElementProblem.hpp" // Constructor - creates the Tpetra objects (maps and vectors) -template typename +template FiniteElementProblem::FiniteElementProblem(int numGlobalElements, - Teuchos::RCP> &comm, - ST s) : - Comm(&comm), - NumGlobalElements(numGlobalElements), - scale(s) -{ - + Teuchos::RCP>& comm, + Scalar s) + : comm(&comm), numGlobalElements(numGlobalElements), scale(s) { // Commonly used variables int i; - MyPID = Comm->getRank(); // Process ID - NumProc = Comm->getSize(); // Total number of processes + 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 tmap_t(NumGlobalElements, 0, *Comm); + standardMap = new map_type(numGlobalElements, 0, *comm); // Get the number of elements owned by this processor - NumMyElements = StandardMap->getLocalNumElements(); + 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 tmap_t(*StandardMap); + 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 OverlapNumMyElements; - int OverlapMinMyGID; - OverlapNumMyElements = NumMyElements + 2; - if ((MyPID == 0) || (MyPID == NumProc - 1)) - OverlapNumMyElements --; - - if (MyPID==0) - OverlapMinMyGID = StandardMap->getMinGlobalIndex(); - else - OverlapMinMyGID = StandardMap->getMinGlobalIndex() - 1; - - int* OverlapMyGlobalElements = new int[OverlapNumMyElements]; + int* overlapMyGlobalElements = new int[overlapNumMyElements]; - for (i = 0; i < OverlapNumMyElements; i ++) - OverlapMyGlobalElements[i] = OverlapMinMyGID + i; + for (i = 0; i < overlapNumMyElements; i++) { + overlapMyGlobalElements[i] = overlapMinMyGID + i; + } - OverlapMap = new tmap_t(-1, OverlapNumMyElements, - OverlapMyGlobalElements, 0, *Comm); + overlapMap = new tmap_t(-1, overlapNumMyElements, overlapMyGlobalElements, 0, *comm); - } // End Overlap map construction ************************************* + } // 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); + 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); @@ -120,41 +117,40 @@ FiniteElementProblem::FiniteElementProblem(int numGlobalElements, } // Matrix and Residual Fills -template typename -bool FiniteElementProblem::evaluate(FillType f, - const tvector_t* soln, - tvector_t* tmp_rhs, - trowmatrix_t* tmp_matrix) -{ +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); + A = dynamic_cast(tmp_matrix); } else if (flag == ALL) { rhs = tmp_rhs; - A = dynamic_cast (tmp_matrix); + 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 - tvector_t u(*OverlapMap); - tvector_t x(*OverlapMap); + vector_type u(*overlapMap); + vector_type x(*overlapMap); // Export Solution to Overlap vector - u.doImport(*soln, *Importer, Tpetra::INSERT); + u.doImport(*soln, *importer, Tpetra::INSERT); // Declare required variables - int i,j,ierr; - int OverlapNumMyElements = OverlapMap->NumMyElements(); + int i, j, ierr; + int overlapNumMyElements = overlapMap->numMyElements(); - int OverlapMinMyGID; - if (MyPID==0) OverlapMinMyGID = StandardMap->getMinGlobalIndex(); - else OverlapMinMyGID = StandardMap->getMinGlobalIndex()-1; + int overlapMinMyGID; + if (myRank == 0) + overlapMinMyGID = standardMap->getMinGlobalIndex(); + else + overlapMinMyGID = standardMap->getMinGlobalIndex() - 1; int row, column; ST jac; @@ -163,152 +159,147 @@ bool FiniteElementProblem::evaluate(FillType f, 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); + 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); + 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++) { - + for (int ne = 0; ne < overlapNumMyElements - 1; ne++) { // Loop Over Gauss Points - for(int gp=0; gp < 2; gp++) { + 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]; + 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", - // MyPID, 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->MyGID(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); + 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 (MyPID==0) { - if ((flag == F_ONLY) || (flag == ALL)) - (*rhs)[0]= (*soln)[0] - leftBC; + 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; + column = 0; + jac = 1.0; A->replaceGlobalValues(0, 1, &jac, &column); - column=1; - jac=0.0; + column = 1; + jac = 0.0; A->replaceGlobalValues(0, 1, &jac, &column); } } - if ( StandardMap->LID(StandardMap->getMaxAllGlobalIndex()) >= 0 ) { - int lastDof = StandardMap->LID(StandardMap->getMaxAllGlobalIndex()); - if ((flag == F_ONLY) || (flag == ALL)) - (*rhs)[lastDof] = (*soln)[lastDof] - rightBC; + 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; + int row = standardMap->getMaxAllGlobalIndex(); + column = row; + jac = 1.0; A->replaceGlobalValues(row, 1, &jac, &column); - column=row-1; - jac=0.0; + column = row - 1; + jac = 0.0; A->replaceGlobalValues(row, 1, &jac, &column); } } // Sync up processors to be safe - Comm->barrier(); + comm->barrier(); A->fillComplete(); return true; } -template typename -Teuchos::RCP FiniteElementProblem::getSolution() -{ +template +Teuchos::RCP FiniteElementProblem::getSolution() { return initialSolution; } -template typename -Teuchos::RCP FiniteElementProblem::getJacobian() -{ +template +Teuchos::RCP FiniteElementProblem::getJacobian() { return A; } -template typename -bool FiniteElementProblem::setParameter(std::string label, ST value) -{ +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; + rightBC = value * scale; else if (label == "Homotopy Continuation Parameter") { // do nothing for now - } - else { + } else { std::cout << "ERROR: FiniteElementProblem::setParameter() - label is invalid " - << "for this problem!" << std::endl; + << "for this problem!" << std::endl; exit(-1); } return true; } -template typename -tcrsgraph_t& FiniteElementProblem::generateGraph(tcrsgraph_t& AAA) -{ +template +tcrsgraph_t& FiniteElementProblem::generateGraph(tcrsgraph_t& AAA) { // Declare required variables - int i,j; + int i, j; int row, column; - int OverlapNumMyElements = OverlapMap->getLocalNumElements(); - int OverlapMinMyGID; - if (MyPID==0) OverlapMinMyGID = StandardMap->getMinGlobalIndex(); - else OverlapMinMyGID = StandardMap->getMinGlobalIndex()-1; + 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++) { - + for (int ne = 0; ne < overlapNumMyElements - 1; ne++) { // Loop over Nodes in Element - for (i=0; i< 2; i++) { - row=OverlapMap->GID(ne+i); + for (i = 0; i < 2; i++) { + row = overlapMap->GID(ne + i); // Loop over Trial Functions - for(j=0;j < 2; j++) { - + 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); + if (standardMap->MyGID(row)) { + column = overlapMap->GID(ne + j); AAA.insertGlobalIndices(row, 1, &column); } } diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.hpp b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.hpp index 70c0378cd437..33c943f8389f 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.hpp +++ b/packages/nox/test/tpetra/LOCA_TestProblems/FiniteElementProblem.hpp @@ -66,41 +66,44 @@ enum FillType {F_ONLY, MATRIX_ONLY, ALL}; // Finite Element Problem Class -template +template class FiniteElementProblem { - using ST = typename Tpetra::Vector::scalar_type; using LO = typename Tpetra::Vector<>::local_ordinal_type; using GO = typename Tpetra::Vector<>::global_ordinal_type; using NO = typename Tpetra::Map<>::node_type; - using tvector_t = typename Tpetra::Vector; - using trowmatrix_t = typename Tpetra::RowMatrix; - using tcrsmatrix_t = typename Tpetra::CrsMatrix; - 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 tvector_t *solnVector, - tvector_t *rhsVector, - trowmatrix_t *matrix, - ST jac_coeff = 1.0, - ST mass_coeff = 0.0) = 0; + 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 tvector_t& getSolution() = 0; + virtual vector_type& getSolution() = 0; // Return a reference to the Tpetra_Vector with the Jacobian // that is generated by the FiniteElementProblem class. - virtual tcrsmatrix_t& getJacobian() = 0; + virtual crs_matrix_type& getJacobian() = 0; // Set a bifurcation parameter in the application physics - virtual bool setParameter(std::string label, ST value) = 0; + virtual bool setParameter(std::string label, scalar_type value) = 0; }; #endif From b03cc067ab84f53a28cb092586f4e067b4ca94c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20Dutheillet-Lamonth=C3=A9zie?= Date: Fri, 29 Sep 2023 16:26:24 +0200 Subject: [PATCH 10/13] #197: add Pitchwork_FiniteElementProblem tpetra version --- .../Pitchfork_FiniteElementProblem.hpp | 153 ++++++++++++++++++ 1 file changed, 153 insertions(+) create mode 100644 packages/nox/test/tpetra/LOCA_TestProblems/Pitchfork_FiniteElementProblem.hpp 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 + + + + From 1200e2f861e57c4571800c7845be40e19fc4aae7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20Dutheillet-Lamonth=C3=A9zie?= Date: Fri, 29 Sep 2023 16:29:11 +0200 Subject: [PATCH 11/13] #197: add Pitchfork definition --- .../Pitchfork_FiniteElementProblem.cpp | 347 ++++++++++++++++++ 1 file changed, 347 insertions(+) create mode 100644 packages/nox/test/tpetra/LOCA_TestProblems/Pitchfork_FiniteElementProblem.cpp 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; +} From 130700fe410f6be2f387bd498e4b567afee9539e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20Dutheillet-Lamonth=C3=A9zie?= Date: Fri, 29 Sep 2023 16:30:01 +0200 Subject: [PATCH 12/13] #197: add LOCALinearConstraint class and new simple test to test classes from the test library --- packages/nox/test/tpetra/CMakeLists.txt | 1 + .../tpetra/LOCA_TestProblems/CMakeLists.txt | 8 +- .../LOCALinearConstraint.cpp | 230 ++++++++++++++++++ .../LOCALinearConstraint.hpp | 199 +++++++++++++++ .../tpetra/LOCA_UnitTests/NGA_DevTests.cpp | 29 +++ 5 files changed, 463 insertions(+), 4 deletions(-) create mode 100644 packages/nox/test/tpetra/LOCA_TestProblems/LOCALinearConstraint.cpp create mode 100644 packages/nox/test/tpetra/LOCA_TestProblems/LOCALinearConstraint.hpp create mode 100644 packages/nox/test/tpetra/LOCA_UnitTests/NGA_DevTests.cpp diff --git a/packages/nox/test/tpetra/CMakeLists.txt b/packages/nox/test/tpetra/CMakeLists.txt index a79074ec3c02..ef652a2ddaab 100644 --- a/packages/nox/test/tpetra/CMakeLists.txt +++ b/packages/nox/test/tpetra/CMakeLists.txt @@ -70,3 +70,4 @@ IF(NOX_ENABLE_ABSTRACT_IMPLEMENTATION_THYRA AND ENDIF() ADD_SUBDIRECTORY(LOCA_TestProblems) +ADD_SUBDIRECTORY(LOCA_UnitTests) diff --git a/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt index 6f2cb03f4396..62035204ecee 100644 --- a/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt +++ b/packages/nox/test/tpetra/LOCA_TestProblems/CMakeLists.txt @@ -9,18 +9,18 @@ APPEND_SET(HEADERS FiniteElementProblem.hpp # Problem_Interface.hpp # Tcubed_FiniteElementProblem.hpp -# Pitchfork_FiniteElementProblem.hpp + Pitchfork_FiniteElementProblem.hpp # --> WIP, TESTING # NormConstraint.hpp -# LOCALinearConstraint.hpp +# LOCALinearConstraint.hpp --> WIP ) APPEND_SET(SOURCES Basis.cpp # Problem_Interface.cpp # Tcubed_FiniteElementProblem.cpp -# Pitchfork_FiniteElementProblem.cpp + Pitchfork_FiniteElementProblem.cpp # --> WIP, TESTING # NormConstraint.cpp -# LOCALinearConstraint.cpp +# LOCALinearConstraint.cpp --> WIP ) # Add this below: 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_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 From ffe144250775fec0dfc53309adb19f9b416d3512 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20Dutheillet-Lamonth=C3=A9zie?= Date: Fri, 29 Sep 2023 16:30:25 +0200 Subject: [PATCH 13/13] #197: update CMakeLists --- .../test/tpetra/LOCA_UnitTests/CMakeLists.txt | 90 +++++++++++++++++++ 1 file changed, 90 insertions(+) create mode 100644 packages/nox/test/tpetra/LOCA_UnitTests/CMakeLists.txt 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()