From 0e1dbb155b3d0e42290f508115c6bbd6def8904e Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Thu, 20 Oct 2022 11:12:05 -0600 Subject: [PATCH] Remove problems and related evaluators (issue #840) --- CMakeLists.txt | 19 - src/CMakeLists.txt | 32 -- src/evaluators/pde/PHAL_CahnHillChemTerm.cpp | 13 - src/evaluators/pde/PHAL_CahnHillChemTerm.hpp | 58 --- .../pde/PHAL_CahnHillChemTerm_Def.hpp | 101 ----- src/evaluators/pde/PHAL_CahnHillRhoResid.cpp | 13 - src/evaluators/pde/PHAL_CahnHillRhoResid.hpp | 68 --- .../pde/PHAL_CahnHillRhoResid_Def.hpp | 125 ------ src/evaluators/pde/PHAL_CahnHillWResid.cpp | 13 - src/evaluators/pde/PHAL_CahnHillWResid.hpp | 60 --- .../pde/PHAL_CahnHillWResid_Def.hpp | 119 ----- src/evaluators/pde/PHAL_ComprNSBodyForce.cpp | 13 - src/evaluators/pde/PHAL_ComprNSBodyForce.hpp | 60 --- .../pde/PHAL_ComprNSBodyForce_Def.hpp | 111 ----- src/evaluators/pde/PHAL_ComprNSResid.cpp | 13 - src/evaluators/pde/PHAL_ComprNSResid.hpp | 76 ---- src/evaluators/pde/PHAL_ComprNSResid_Def.hpp | 202 --------- src/evaluators/pde/PHAL_ComprNSViscosity.cpp | 13 - src/evaluators/pde/PHAL_ComprNSViscosity.hpp | 75 ---- .../pde/PHAL_ComprNSViscosity_Def.hpp | 178 -------- .../pde/PHAL_LinComprNSBodyForce.cpp | 13 - .../pde/PHAL_LinComprNSBodyForce.hpp | 60 --- .../pde/PHAL_LinComprNSBodyForce_Def.hpp | 143 ------ src/evaluators/pde/PHAL_LinComprNSResid.cpp | 13 - src/evaluators/pde/PHAL_LinComprNSResid.hpp | 72 --- .../pde/PHAL_LinComprNSResid_Def.hpp | 422 ------------------ src/evaluators/pde/PNP_ConcentrationResid.cpp | 13 - src/evaluators/pde/PNP_ConcentrationResid.hpp | 61 --- .../pde/PNP_ConcentrationResid_Def.hpp | 104 ----- src/evaluators/pde/PNP_PotentialResid.cpp | 13 - src/evaluators/pde/PNP_PotentialResid.hpp | 58 --- src/evaluators/pde/PNP_PotentialResid_Def.hpp | 92 ---- src/problems/Albany_CahnHillProblem.cpp | 105 ----- src/problems/Albany_CahnHillProblem.hpp | 318 ------------- src/problems/Albany_ComprNSProblem.cpp | 200 --------- src/problems/Albany_ComprNSProblem.hpp | 314 ------------- src/problems/Albany_DemoProblemFactory.cpp | 35 +- src/problems/Albany_LinComprNSProblem.cpp | 99 ---- src/problems/Albany_LinComprNSProblem.hpp | 272 ----------- src/problems/Albany_PNPProblem.cpp | 178 -------- src/problems/Albany_PNPProblem.hpp | 304 ------------- src/problems/Albany_ProblemUtils.cpp | 8 - 42 files changed, 2 insertions(+), 4257 deletions(-) delete mode 100644 src/evaluators/pde/PHAL_CahnHillChemTerm.cpp delete mode 100644 src/evaluators/pde/PHAL_CahnHillChemTerm.hpp delete mode 100644 src/evaluators/pde/PHAL_CahnHillChemTerm_Def.hpp delete mode 100644 src/evaluators/pde/PHAL_CahnHillRhoResid.cpp delete mode 100644 src/evaluators/pde/PHAL_CahnHillRhoResid.hpp delete mode 100644 src/evaluators/pde/PHAL_CahnHillRhoResid_Def.hpp delete mode 100644 src/evaluators/pde/PHAL_CahnHillWResid.cpp delete mode 100644 src/evaluators/pde/PHAL_CahnHillWResid.hpp delete mode 100644 src/evaluators/pde/PHAL_CahnHillWResid_Def.hpp delete mode 100644 src/evaluators/pde/PHAL_ComprNSBodyForce.cpp delete mode 100644 src/evaluators/pde/PHAL_ComprNSBodyForce.hpp delete mode 100644 src/evaluators/pde/PHAL_ComprNSBodyForce_Def.hpp delete mode 100644 src/evaluators/pde/PHAL_ComprNSResid.cpp delete mode 100644 src/evaluators/pde/PHAL_ComprNSResid.hpp delete mode 100644 src/evaluators/pde/PHAL_ComprNSResid_Def.hpp delete mode 100644 src/evaluators/pde/PHAL_ComprNSViscosity.cpp delete mode 100644 src/evaluators/pde/PHAL_ComprNSViscosity.hpp delete mode 100644 src/evaluators/pde/PHAL_ComprNSViscosity_Def.hpp delete mode 100644 src/evaluators/pde/PHAL_LinComprNSBodyForce.cpp delete mode 100644 src/evaluators/pde/PHAL_LinComprNSBodyForce.hpp delete mode 100644 src/evaluators/pde/PHAL_LinComprNSBodyForce_Def.hpp delete mode 100644 src/evaluators/pde/PHAL_LinComprNSResid.cpp delete mode 100644 src/evaluators/pde/PHAL_LinComprNSResid.hpp delete mode 100644 src/evaluators/pde/PHAL_LinComprNSResid_Def.hpp delete mode 100644 src/evaluators/pde/PNP_ConcentrationResid.cpp delete mode 100644 src/evaluators/pde/PNP_ConcentrationResid.hpp delete mode 100644 src/evaluators/pde/PNP_ConcentrationResid_Def.hpp delete mode 100644 src/evaluators/pde/PNP_PotentialResid.cpp delete mode 100644 src/evaluators/pde/PNP_PotentialResid.hpp delete mode 100644 src/evaluators/pde/PNP_PotentialResid_Def.hpp delete mode 100644 src/problems/Albany_CahnHillProblem.cpp delete mode 100644 src/problems/Albany_CahnHillProblem.hpp delete mode 100644 src/problems/Albany_ComprNSProblem.cpp delete mode 100644 src/problems/Albany_ComprNSProblem.hpp delete mode 100644 src/problems/Albany_LinComprNSProblem.cpp delete mode 100644 src/problems/Albany_LinComprNSProblem.hpp delete mode 100644 src/problems/Albany_PNPProblem.cpp delete mode 100644 src/problems/Albany_PNPProblem.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 1716f06035..e95bf04dd5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -571,16 +571,6 @@ ELSE() SET(ALBANY_MESH_DEPENDS_ON_PARAMETERS FALSE) ENDIF() -# set optional dependency of mesh on solution, defaults to Disabled -OPTION(ENABLE_MESH_DEPENDS_ON_SOLUTION "Flag to turn on dependency of mesh on solution" OFF) -IF (ENABLE_MESH_DEPENDS_ON_SOLUTION) - MESSAGE("-- MESH_DEPENDS_ON_SOLUTION Enabled.") - SET(ALBANY_MESH_DEPENDS_ON_SOLUTION TRUE) -ELSE() - MESSAGE("-- MESH_DEPENDS_ON_SOLUTION NOT Enabled.") - SET(ALBANY_MESH_DEPENDS_ON_SOLUTION FALSE) -ENDIF() - MESSAGE("\nAlbany physics packages:\n") # set optional dependency on demoPDEs, defaults to Enabled @@ -734,15 +724,6 @@ ELSE() SET(ALBANY_FADTYPE_NOTEQUAL_TANFADTYPE TRUE) MESSAGE("-- FAD_TYPE is not TAN_FAD_TYPE") ENDIF() - -LIST(FIND Trilinos_PACKAGE_LIST Pamgen PAMGEN_List_ID) -IF (NOT PAMGEN_List_ID GREATER -1) - MESSAGE("-- Pamgen is Enabled. Building Pamgen tests") - set(ALBANY_PAMGEN TRUE) -ELSE() - MESSAGE("-- Pamgen is NOT Enabled. Not building Pamgen tests") - set(ALBANY_PAMGEN FALSE) -ENDIF() # Disable the RTC capability if Trilinos is not built with Pamgen LIST(FIND Trilinos_PACKAGE_LIST Pamgen PAMGEN_List_ID) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c91a4687ed..30110e42a7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -615,29 +615,17 @@ if(ALBANY_DEMO_PDES) set(PDE_SRCS problems/Albany_AdvDiffProblem.cpp problems/Albany_ReactDiffSystem.cpp - problems/Albany_CahnHillProblem.cpp - problems/Albany_ComprNSProblem.cpp problems/Albany_DemoProblemFactory.cpp problems/Albany_Helmholtz2DProblem.cpp - problems/Albany_LinComprNSProblem.cpp problems/Albany_ODEProblem.cpp problems/Albany_NavierStokes.cpp - problems/Albany_PNPProblem.cpp problems/Albany_ThermoElectrostaticsProblem.cpp problems/Albany_ThermalProblem.cpp problems/Albany_AdvectionProblem.cpp evaluators/pde/PHAL_AdvDiffResid.cpp evaluators/pde/PHAL_ReactDiffSystemResid.cpp - evaluators/pde/PHAL_CahnHillChemTerm.cpp - evaluators/pde/PHAL_CahnHillRhoResid.cpp - evaluators/pde/PHAL_CahnHillWResid.cpp - evaluators/pde/PHAL_ComprNSBodyForce.cpp - evaluators/pde/PHAL_ComprNSResid.cpp - evaluators/pde/PHAL_ComprNSViscosity.cpp evaluators/pde/PHAL_HelmholtzResid.cpp evaluators/pde/PHAL_JouleHeating.cpp - evaluators/pde/PHAL_LinComprNSBodyForce.cpp - evaluators/pde/PHAL_LinComprNSResid.cpp evaluators/pde/PHAL_NSBodyForce.cpp evaluators/pde/PHAL_NSContinuityResid.cpp evaluators/pde/PHAL_NSContravarientMetricTensor.cpp @@ -655,20 +643,14 @@ if(ALBANY_DEMO_PDES) evaluators/pde/PHAL_PoissonResid.cpp evaluators/pde/PHAL_Permittivity.cpp evaluators/pde/PHAL_TEProp.cpp - evaluators/pde/PNP_ConcentrationResid.cpp - evaluators/pde/PNP_PotentialResid.cpp ) set(PDE_HDRS problems/Albany_AdvDiffProblem.hpp problems/Albany_ReactDiffSystem.hpp - problems/Albany_CahnHillProblem.hpp - problems/Albany_ComprNSProblem.hpp problems/Albany_Helmholtz2DProblem.hpp - problems/Albany_LinComprNSProblem.hpp problems/Albany_NavierStokes.hpp problems/Albany_ODEProblem.hpp - problems/Albany_PNPProblem.hpp problems/Albany_ThermoElectrostaticsProblem.hpp problems/Albany_ThermalProblem.hpp problems/Albany_AdvectionProblem.hpp @@ -676,12 +658,6 @@ if(ALBANY_DEMO_PDES) evaluators/pde/PHAL_AdvDiffResid_Def.hpp evaluators/pde/PHAL_ReactDiffSystemResid.hpp evaluators/pde/PHAL_ReactDiffSystemResid_Def.hpp - evaluators/pde/PHAL_CahnHillRhoResid.hpp - evaluators/pde/PHAL_CahnHillRhoResid_Def.hpp - evaluators/pde/PHAL_CahnHillWResid.hpp - evaluators/pde/PHAL_CahnHillWResid_Def.hpp - evaluators/pde/PHAL_CahnHillChemTerm.hpp - evaluators/pde/PHAL_CahnHillChemTerm_Def.hpp evaluators/pde/PHAL_ComprNSBodyForce.hpp evaluators/pde/PHAL_ComprNSBodyForce_Def.hpp evaluators/pde/PHAL_ComprNSResid.hpp @@ -692,10 +668,6 @@ if(ALBANY_DEMO_PDES) evaluators/pde/PHAL_HelmholtzResid_Def.hpp evaluators/pde/PHAL_JouleHeating.hpp evaluators/pde/PHAL_JouleHeating_Def.hpp - evaluators/pde/PHAL_LinComprNSBodyForce.hpp - evaluators/pde/PHAL_LinComprNSBodyForce_Def.hpp - evaluators/pde/PHAL_LinComprNSResid.hpp - evaluators/pde/PHAL_LinComprNSResid_Def.hpp evaluators/pde/PHAL_NSContinuityResid.hpp evaluators/pde/PHAL_NSContinuityResid_Def.hpp evaluators/pde/PHAL_NSBodyForce.hpp @@ -726,10 +698,6 @@ if(ALBANY_DEMO_PDES) evaluators/pde/PHAL_Permittivity_Def.hpp evaluators/pde/PHAL_TEProp.hpp evaluators/pde/PHAL_TEProp_Def.hpp - evaluators/pde/PNP_ConcentrationResid.hpp - evaluators/pde/PNP_ConcentrationResid_Def.hpp - evaluators/pde/PNP_PotentialResid.hpp - evaluators/pde/PNP_PotentialResid_Def.hpp evaluators/pde/PHAL_ThermalResid.hpp evaluators/pde/PHAL_ThermalResid_Def.hpp evaluators/pde/PHAL_AdvectionResid.hpp diff --git a/src/evaluators/pde/PHAL_CahnHillChemTerm.cpp b/src/evaluators/pde/PHAL_CahnHillChemTerm.cpp deleted file mode 100644 index 65c800a047..0000000000 --- a/src/evaluators/pde/PHAL_CahnHillChemTerm.cpp +++ /dev/null @@ -1,13 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "PHAL_AlbanyTraits.hpp" - -#include "PHAL_CahnHillChemTerm.hpp" -#include "PHAL_CahnHillChemTerm_Def.hpp" - -PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::CahnHillChemTerm) - diff --git a/src/evaluators/pde/PHAL_CahnHillChemTerm.hpp b/src/evaluators/pde/PHAL_CahnHillChemTerm.hpp deleted file mode 100644 index dcedc4d7a1..0000000000 --- a/src/evaluators/pde/PHAL_CahnHillChemTerm.hpp +++ /dev/null @@ -1,58 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef PHAL_CAHNHILLCHEMTERM_HPP -#define PHAL_CAHNHILLCHEMTERM_HPP - -#include "Phalanx_config.hpp" -#include "Phalanx_Evaluator_WithBaseImpl.hpp" -#include "Phalanx_Evaluator_Derived.hpp" -#include "Phalanx_MDField.hpp" - -/** \brief Finite Element Interpolation Evaluator - - This evaluator interpolates nodal DOF values to quad points. - -*/ -namespace PHAL { - -template -class CahnHillChemTerm : public PHX::EvaluatorWithBaseImpl, - public PHX::EvaluatorDerived { - -public: - - typedef typename EvalT::ScalarT ScalarT; - - CahnHillChemTerm(const Teuchos::ParameterList& p); - - void postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& vm); - - void evaluateFields(typename Traits::EvalData d); - - ScalarT& getValue(const std::string &n); - - -private: - - typedef typename EvalT::MeshScalarT MeshScalarT; - - // Input: - PHX::MDField rho; - PHX::MDField w; - - // Output: - PHX::MDField chemTerm; - - unsigned int numQPs, numDims, numNodes; - - ScalarT b; - -}; -} - -#endif diff --git a/src/evaluators/pde/PHAL_CahnHillChemTerm_Def.hpp b/src/evaluators/pde/PHAL_CahnHillChemTerm_Def.hpp deleted file mode 100644 index a16dac9bff..0000000000 --- a/src/evaluators/pde/PHAL_CahnHillChemTerm_Def.hpp +++ /dev/null @@ -1,101 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Teuchos_TestForException.hpp" -#include "Phalanx_DataLayout.hpp" - -#include "Intrepid2_FunctionSpaceTools.hpp" - - -template -inline ScalarT Sqr (const ScalarT& num) { - return num * num; -} - -namespace PHAL { - - -//********************************************************************** -template -CahnHillChemTerm:: -CahnHillChemTerm(const Teuchos::ParameterList& p) : - rho (p.get ("Rho QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - w (p.get ("W QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - chemTerm (p.get ("Chemical Energy Term"), - p.get >("QP Scalar Data Layout") ) - -{ - - b = p.get("b Value"); - - Teuchos::RCP vector_dl = - p.get< Teuchos::RCP >("QP Vector Data Layout"); - std::vector dims; - vector_dl->dimensions(dims); - numQPs = dims[1]; - numDims = dims[2]; - - this->addDependentField(rho.fieldTag()); - this->addDependentField(w.fieldTag()); - - this->addEvaluatedField(chemTerm); - - this->setName("CahnHillChemTerm" ); - -} - -//********************************************************************** -template -void CahnHillChemTerm:: -postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& fm) -{ - this->utils.setFieldData(rho,fm); - this->utils.setFieldData(w,fm); - - this->utils.setFieldData(chemTerm,fm); -} - -//********************************************************************** -template -void CahnHillChemTerm:: -evaluateFields(typename Traits::EvalData workset) -{ - -// Equations 1.1 and 2.2 in Garcke, Rumpf, and Weikard -// psi(rho) = 0.25 * (rho^2 - b^2)^2 - - for (std::size_t cell=0; cell < workset.numCells; ++cell) - for (std::size_t qp=0; qp < numQPs; ++qp) - - // chemTerm(cell, qp) = 0.25 * Sqr(Sqr(rho(cell, qp)) - Sqr(b)) - w(cell, qp); - chemTerm(cell, qp) = ( Sqr(rho(cell, qp)) - Sqr(b) ) * rho(cell, qp) - w(cell, qp); - -} - -template -typename CahnHillChemTerm::ScalarT& -CahnHillChemTerm::getValue(const std::string &n) { - - if (n == "b") - - return b; - - else - { - TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, std::endl << - "Error! Logic error in getting parameter " << n << - " in CahnHillChemTerm::getValue()" << std::endl); - return b; - } - -} - - -} - diff --git a/src/evaluators/pde/PHAL_CahnHillRhoResid.cpp b/src/evaluators/pde/PHAL_CahnHillRhoResid.cpp deleted file mode 100644 index f0d1901f7c..0000000000 --- a/src/evaluators/pde/PHAL_CahnHillRhoResid.cpp +++ /dev/null @@ -1,13 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "PHAL_AlbanyTraits.hpp" - -#include "PHAL_CahnHillRhoResid.hpp" -#include "PHAL_CahnHillRhoResid_Def.hpp" - -PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::CahnHillRhoResid) - diff --git a/src/evaluators/pde/PHAL_CahnHillRhoResid.hpp b/src/evaluators/pde/PHAL_CahnHillRhoResid.hpp deleted file mode 100644 index 1a736e4733..0000000000 --- a/src/evaluators/pde/PHAL_CahnHillRhoResid.hpp +++ /dev/null @@ -1,68 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef PHAL_CAHNHILLRHORESID_HPP -#define PHAL_CAHNHILLRHORESID_HPP - -#include "Phalanx_config.hpp" -#include "Phalanx_Evaluator_WithBaseImpl.hpp" -#include "Phalanx_Evaluator_Derived.hpp" -#include "Phalanx_MDField.hpp" - -namespace PHAL { - -/** \brief Finite Element Interpolation Evaluator - - This evaluator interpolates nodal DOF values to quad points. - -*/ - -template -class CahnHillRhoResid : public PHX::EvaluatorWithBaseImpl, - public PHX::EvaluatorDerived { - -public: - - typedef typename EvalT::ScalarT ScalarT; - - CahnHillRhoResid(const Teuchos::ParameterList& p); - - void postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& vm); - - void evaluateFields(typename Traits::EvalData d); - - ScalarT& getValue(const std::string &n); - -private: - - typedef typename EvalT::MeshScalarT MeshScalarT; - - // Input: - PHX::MDField wBF; - PHX::MDField wGradBF; - PHX::MDField rhoGrad; - PHX::MDField chemTerm; - PHX::MDField noiseTerm; - - - // Output: - PHX::MDField rhoResidual; - - Kokkos::DynRankView gamma_term; - - unsigned int numQPs, numDims, numNodes, worksetSize; - - ScalarT gamma; - - // Langevin noise present - bool haveNoise; - - -}; -} - -#endif diff --git a/src/evaluators/pde/PHAL_CahnHillRhoResid_Def.hpp b/src/evaluators/pde/PHAL_CahnHillRhoResid_Def.hpp deleted file mode 100644 index 9c485e8414..0000000000 --- a/src/evaluators/pde/PHAL_CahnHillRhoResid_Def.hpp +++ /dev/null @@ -1,125 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Teuchos_TestForException.hpp" -#include "Phalanx_DataLayout.hpp" - -#include "Intrepid2_FunctionSpaceTools.hpp" - -namespace PHAL { - -//********************************************************************** -template -CahnHillRhoResid:: -CahnHillRhoResid(const Teuchos::ParameterList& p) : - wBF (p.get ("Weighted BF Name"), - p.get >("Node QP Scalar Data Layout") ), - wGradBF (p.get ("Weighted Gradient BF Name"), - p.get >("Node QP Vector Data Layout") ), - rhoGrad (p.get ("Gradient QP Variable Name"), - p.get >("QP Vector Data Layout") ), - chemTerm (p.get ("Chemical Energy Term"), - p.get >("QP Scalar Data Layout") ), - rhoResidual (p.get ("Residual Name"), - p.get >("Node Scalar Data Layout") ) -{ - - haveNoise = p.get("Have Noise"); - - this->addDependentField(wBF.fieldTag()); - this->addDependentField(rhoGrad.fieldTag()); - this->addDependentField(wGradBF.fieldTag()); - this->addDependentField(chemTerm.fieldTag()); - this->addEvaluatedField(rhoResidual); - - if(haveNoise){ - noiseTerm = decltype(noiseTerm)(p.get("Langevin Noise Term"), - p.get >("QP Scalar Data Layout") ), - this->addDependentField(noiseTerm.fieldTag()); - } - - gamma = p.get("gamma Value"); - - Teuchos::RCP vector_dl = - p.get< Teuchos::RCP >("Node QP Vector Data Layout"); - std::vector dims; - vector_dl->dimensions(dims); - worksetSize = dims[0]; - numNodes = dims[1]; - numQPs = dims[2]; - numDims = dims[3]; - - this->setName("CahnHillRhoResid" ); - -} - -//********************************************************************** -template -void CahnHillRhoResid:: -postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& fm) -{ - this->utils.setFieldData(wBF,fm); - this->utils.setFieldData(rhoGrad,fm); - this->utils.setFieldData(wGradBF,fm); - this->utils.setFieldData(chemTerm,fm); - if(haveNoise) - this->utils.setFieldData(noiseTerm,fm); - - this->utils.setFieldData(rhoResidual,fm); - - gamma_term = Kokkos::createDynRankView(rhoGrad.get_view(), "XXX", worksetSize, numQPs, numDims); - -} - -//********************************************************************** -template -void CahnHillRhoResid:: -evaluateFields(typename Traits::EvalData workset) -{ - -// Form Equation 2.2 - - typedef Intrepid2::FunctionSpaceTools FST; - - for (std::size_t cell=0; cell < workset.numCells; ++cell) - for (std::size_t qp=0; qp < numQPs; ++qp) - for (std::size_t i=0; i < numDims; ++i) - - gamma_term(cell, qp, i) = rhoGrad(cell,qp,i) * gamma; - - FST::integrate(rhoResidual.get_view(), gamma_term, wGradBF.get_view(), false); // "false" overwrites - - FST::integrate(rhoResidual.get_view(), chemTerm.get_view(), wBF.get_view(), true); // "true" sums into - - if(haveNoise) - - FST::integrate(rhoResidual.get_view(), noiseTerm.get_view(), wBF.get_view(), true); // "true" sums into - - -} - -template -typename CahnHillRhoResid::ScalarT& -CahnHillRhoResid::getValue(const std::string &n) { - - if (n == "gamma") - - return gamma; - - else - { - TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, std::endl << - "Error! Logic error in getting parameter " << n << - " in CahnHillRhoResid::getValue()" << std::endl); - return gamma; - } - -} - -//********************************************************************** -} - diff --git a/src/evaluators/pde/PHAL_CahnHillWResid.cpp b/src/evaluators/pde/PHAL_CahnHillWResid.cpp deleted file mode 100644 index c03443963c..0000000000 --- a/src/evaluators/pde/PHAL_CahnHillWResid.cpp +++ /dev/null @@ -1,13 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "PHAL_AlbanyTraits.hpp" - -#include "PHAL_CahnHillWResid.hpp" -#include "PHAL_CahnHillWResid_Def.hpp" - -PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::CahnHillWResid) - diff --git a/src/evaluators/pde/PHAL_CahnHillWResid.hpp b/src/evaluators/pde/PHAL_CahnHillWResid.hpp deleted file mode 100644 index 97a14e8e28..0000000000 --- a/src/evaluators/pde/PHAL_CahnHillWResid.hpp +++ /dev/null @@ -1,60 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef PHAL_CAHNHILLWRESID_HPP -#define PHAL_CAHNHILLWRESID_HPP - -#include "Phalanx_config.hpp" -#include "Phalanx_Evaluator_WithBaseImpl.hpp" -#include "Phalanx_Evaluator_Derived.hpp" -#include "Phalanx_MDField.hpp" - -namespace PHAL { - -/** \brief Finite Element Interpolation Evaluator - - This evaluator interpolates nodal DOF values to quad points. - -*/ - -template -class CahnHillWResid : public PHX::EvaluatorWithBaseImpl, - public PHX::EvaluatorDerived { - -public: - - CahnHillWResid(const Teuchos::ParameterList& p); - - void postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& vm); - - void evaluateFields(typename Traits::EvalData d); - -private: - - typedef typename EvalT::ScalarT ScalarT; - typedef typename EvalT::MeshScalarT MeshScalarT; - - // Input: - PHX::MDField wBF; - PHX::MDField BF; - PHX::MDField rhoDot; - PHX::MDField rhoDotNode; - PHX::MDField wGradBF; - PHX::MDField wGrad; - - // Output: - PHX::MDField wResidual; - - unsigned int numQPs, numDims, numNodes, worksetSize; - - // lump mass matrix - bool lump; - -}; -} - -#endif diff --git a/src/evaluators/pde/PHAL_CahnHillWResid_Def.hpp b/src/evaluators/pde/PHAL_CahnHillWResid_Def.hpp deleted file mode 100644 index 3338c9bd6c..0000000000 --- a/src/evaluators/pde/PHAL_CahnHillWResid_Def.hpp +++ /dev/null @@ -1,119 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Teuchos_TestForException.hpp" -#include "Phalanx_DataLayout.hpp" - -#include "Intrepid2_FunctionSpaceTools.hpp" - -namespace PHAL { - -//********************************************************************** -template -CahnHillWResid:: -CahnHillWResid(const Teuchos::ParameterList& p) : - wBF (p.get ("Weighted BF Name"), - p.get >("Node QP Scalar Data Layout") ), - BF (p.get ("BF Name"), - p.get >("Node QP Scalar Data Layout") ), - rhoDot (p.get ("Rho QP Time Derivative Variable Name"), - p.get >("QP Scalar Data Layout") ), - rhoDotNode (p.get ("Rho QP Time Derivative Variable Name"), - p.get >("Node Scalar Data Layout") ), - wGradBF (p.get ("Weighted Gradient BF Name"), - p.get >("Node QP Vector Data Layout") ), - wGrad (p.get ("Gradient QP Variable Name"), - p.get >("QP Vector Data Layout") ), - wResidual (p.get ("Residual Name"), - p.get >("Node Scalar Data Layout") ) -{ - - lump = p.get("Lump Mass"); - - this->addDependentField(wBF.fieldTag()); - this->addDependentField(BF.fieldTag()); - this->addDependentField(rhoDot.fieldTag()); - this->addDependentField(rhoDotNode.fieldTag()); - this->addDependentField(wGrad.fieldTag()); - this->addDependentField(wGradBF.fieldTag()); - this->addEvaluatedField(wResidual); - - Teuchos::RCP vector_dl = - p.get< Teuchos::RCP >("Node QP Vector Data Layout"); - std::vector dims; - vector_dl->dimensions(dims); - worksetSize = dims[0]; - numNodes = dims[1]; - numQPs = dims[2]; - numDims = dims[3]; - - this->setName("CahnHillWResid" ); - -} - -//********************************************************************** -template -void CahnHillWResid:: -postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& fm) -{ - this->utils.setFieldData(wBF,fm); - this->utils.setFieldData(BF,fm); - this->utils.setFieldData(wGrad,fm); - this->utils.setFieldData(wGradBF,fm); - this->utils.setFieldData(rhoDot,fm); - this->utils.setFieldData(rhoDotNode,fm); - - this->utils.setFieldData(wResidual,fm); -} - -template -void CahnHillWResid:: -evaluateFields(typename Traits::EvalData workset) -{ - typedef Intrepid2::FunctionSpaceTools FST; - - FST::integrate(wResidual.get_view(), wGrad.get_view(), wGradBF.get_view(), false); // "false" overwrites - - if(!lump){ - // Consistent mass matrix, the Intrepid2 way - FST::integrate(wResidual.get_view(), rhoDot.get_view(), wBF.get_view(), true); // "true" sums into - - // Consistent mass matrix, done manually -/* - for (std::size_t cell=0; cell < workset.numCells; ++cell) - for (std::size_t node=0; node < numNodes; ++node) - for (std::size_t qp=0; qp < numQPs; ++qp) - - wResidual(cell, node) += rhoDot(cell, qp) * wBF(cell, node, qp); -*/ - } - else { - - RealType diag; - - // Lumped mass matrix - for (std::size_t cell=0; cell < workset.numCells; ++cell) - for (std::size_t qp=0; qp < numQPs; ++qp) { - - diag = 0; - for (std::size_t node=0; node < numNodes; ++node) - - diag += BF(cell, node, qp); // lump all the row onto the diagonal - - for (std::size_t node=0; node < numNodes; ++node) - - wResidual(cell, node) += diag * rhoDotNode(cell, node) * wBF(cell, node, qp); - - } - - } - -} - -//********************************************************************** -} - diff --git a/src/evaluators/pde/PHAL_ComprNSBodyForce.cpp b/src/evaluators/pde/PHAL_ComprNSBodyForce.cpp deleted file mode 100644 index 6ca36a0270..0000000000 --- a/src/evaluators/pde/PHAL_ComprNSBodyForce.cpp +++ /dev/null @@ -1,13 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "PHAL_AlbanyTraits.hpp" - -#include "PHAL_ComprNSBodyForce.hpp" -#include "PHAL_ComprNSBodyForce_Def.hpp" - -PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::ComprNSBodyForce) - diff --git a/src/evaluators/pde/PHAL_ComprNSBodyForce.hpp b/src/evaluators/pde/PHAL_ComprNSBodyForce.hpp deleted file mode 100644 index a65e9d487f..0000000000 --- a/src/evaluators/pde/PHAL_ComprNSBodyForce.hpp +++ /dev/null @@ -1,60 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef PHAL_COMPRNSBODYFORCE_HPP -#define PHAL_COMPRNSBODYFORCE_HPP - -#include "Phalanx_config.hpp" -#include "Phalanx_Evaluator_WithBaseImpl.hpp" -#include "Phalanx_Evaluator_Derived.hpp" -#include "Phalanx_MDField.hpp" - -namespace PHAL { -/** \brief Finite Element Interpolation Evaluator - - This evaluator interpolates nodal DOF values to quad points. - -*/ - -template -class ComprNSBodyForce : public PHX::EvaluatorWithBaseImpl, - public PHX::EvaluatorDerived { - -public: - - typedef typename EvalT::ScalarT ScalarT; - - ComprNSBodyForce(const Teuchos::ParameterList& p); - - void postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& vm); - - void evaluateFields(typename Traits::EvalData d); - - -private: - - typedef typename EvalT::MeshScalarT MeshScalarT; - - // Input: - PHX::MDField coordVec; - Teuchos::Array gravity; - - // Output: - PHX::MDField force; - - //Force types - enum BFTYPE {NONE, TAYLOR_GREEN_VORTEX}; - BFTYPE bf_type; - - std::size_t numQPs; - std::size_t numDims; - std::size_t vecDim; - -}; -} - -#endif diff --git a/src/evaluators/pde/PHAL_ComprNSBodyForce_Def.hpp b/src/evaluators/pde/PHAL_ComprNSBodyForce_Def.hpp deleted file mode 100644 index d782e747cf..0000000000 --- a/src/evaluators/pde/PHAL_ComprNSBodyForce_Def.hpp +++ /dev/null @@ -1,111 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Teuchos_TestForException.hpp" -#include "Phalanx_DataLayout.hpp" -#include "Sacado.hpp" - -#include "Intrepid2_FunctionSpaceTools.hpp" - -namespace PHAL { -const double pi = 3.1415926535897932385; -//********************************************************************** - -template -ComprNSBodyForce:: -ComprNSBodyForce(const Teuchos::ParameterList& p) : - force(p.get("Body Force Name"), - p.get >("QP Vector Data Layout") ) -{ - std::cout << "Compr NS body force constructor!" << std::endl; - Teuchos::ParameterList* bf_list = - p.get("Parameter List"); - - std::string type = bf_list->get("Type", "None"); - if (type == "None") { - bf_type = NONE; - } - else if (type == "Taylor-Green Vortex") { - std::cout << "Taylor-Green Vortex source" << std::endl; - bf_type = TAYLOR_GREEN_VORTEX; - coordVec = decltype(coordVec)( - p.get("Coordinate Vector Name"), - p.get >("QP Gradient Data Layout") ); - this->addDependentField(coordVec.fieldTag()); - } - - this->addEvaluatedField(force); - - Teuchos::RCP gradient_dl = - p.get< Teuchos::RCP >("QP Gradient Data Layout"); - std::vector dims; - gradient_dl->dimensions(dims); - numQPs = dims[1]; - numDims = dims[2]; - Teuchos::RCP vector_dl = - p.get< Teuchos::RCP >("QP Vector Data Layout"); - vector_dl->dimensions(dims); - vecDim = dims[2]; - -std::cout << " in Compr NS Stokes source! " << std::endl; -std::cout << " vecDim = " << vecDim << std::endl; -std::cout << " numDims = " << numDims << std::endl; -std::cout << " numQPs = " << numQPs << std::endl; - - - this->setName("ComprNSBodyForce" ); -} - -//********************************************************************** -template -void ComprNSBodyForce:: -postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& fm) -{ - if (bf_type == TAYLOR_GREEN_VORTEX) { - this->utils.setFieldData(coordVec,fm); - } - this->utils.setFieldData(force,fm); -} - -//********************************************************************** -template -void ComprNSBodyForce:: -evaluateFields(typename Traits::EvalData workset) -{ - if (bf_type == NONE) { - for (std::size_t cell=0; cell < workset.numCells; ++cell) - for (std::size_t qp=0; qp < numQPs; ++qp) - for (std::size_t i=0; i < vecDim; ++i) - force(cell,qp,i) = 0.0; - } - else if (bf_type == TAYLOR_GREEN_VORTEX) { //source term for MMS Taylor-Vortex-like problem - const RealType time = workset.current_time; //time - const double Re = 1.0; - const double Pr = 0.72; - const double gamma_gas = 1.4; - const double kappa = 1.0; - const double mu = 1.0/Re; - const double Rgas = 0.714285733; //non-dimensional gas constant - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t qp=0; qp < numQPs; ++qp) { - MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0); - MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1); - force(cell,qp,0) = 0.0; - force(cell,qp,1) = 2.0*exp(-2.0*time)*cos(x2pi)*(-sin(y2pi) + exp(-2.0*time)*sin(x2pi)*pi + 4.0*mu*pi*pi*sin(y2pi)) + 2.0*Rgas*pi*exp(-4.0*time)*sin(x2pi); - force(cell,qp,2) = 2.0*exp(-2.0*time)*cos(y2pi)*(sin(x2pi) + exp(-2.0*time)*sin(y2pi)*pi - 4.0*mu*pi*pi*sin(x2pi)) + 2.0*Rgas*pi*exp(-4.0*time)*sin(y2pi); - force(cell,qp,3) = -2.0*exp(-4.0*time)*(-2.0*cos(x2pi) - 2.0*cos(y2pi) + exp(-2.0*time)*cos(x2pi)*sin(y2pi)*sin(x2pi)*pi - - exp(-2.0*time)*sin(x2pi)*cos(y2pi)*sin(y2pi)*pi) - + (gamma_gas-1.0)/Rgas*4.0*mu*exp(-2.0*time)*sin(x2pi)*sin(y2pi)*pi*2.0*exp(-2.0*time)*sin(x2pi)*pi*sin(y2pi) - + (gamma_gas-1.0)/Rgas*4.0*mu*exp(-2.0*time)*sin(x2pi)*pi*sin(y2pi)*2.0*exp(-2.0*time)*sin(x2pi)*pi*sin(y2pi) - - gamma_gas*kappa/(Pr*Re)*4.0*exp(-4.0*time)*pi*pi*(cos(x2pi) + cos(y2pi)); - } - } - } -} - -} - diff --git a/src/evaluators/pde/PHAL_ComprNSResid.cpp b/src/evaluators/pde/PHAL_ComprNSResid.cpp deleted file mode 100644 index 12ccb3afa4..0000000000 --- a/src/evaluators/pde/PHAL_ComprNSResid.cpp +++ /dev/null @@ -1,13 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "PHAL_AlbanyTraits.hpp" - -#include "PHAL_ComprNSResid.hpp" -#include "PHAL_ComprNSResid_Def.hpp" - -PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::ComprNSResid) - diff --git a/src/evaluators/pde/PHAL_ComprNSResid.hpp b/src/evaluators/pde/PHAL_ComprNSResid.hpp deleted file mode 100644 index 337b0c6945..0000000000 --- a/src/evaluators/pde/PHAL_ComprNSResid.hpp +++ /dev/null @@ -1,76 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef PHAL_COMPRNSRESID_HPP -#define PHAL_COMPRNSRESID_HPP - -#include "Phalanx_config.hpp" -#include "Phalanx_Evaluator_WithBaseImpl.hpp" -#include "Phalanx_Evaluator_Derived.hpp" -#include "Phalanx_MDField.hpp" - -namespace PHAL { -/** \brief Finite Element Interpolation Evaluator - - This evaluator interpolates nodal DOF values to quad points. - -*/ - -template -class ComprNSResid : public PHX::EvaluatorWithBaseImpl, - public PHX::EvaluatorDerived { - -public: - - ComprNSResid(const Teuchos::ParameterList& p); - - void postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& vm); - - void evaluateFields(typename Traits::EvalData d); - -private: - - typedef typename EvalT::ScalarT ScalarT; - typedef typename EvalT::MeshScalarT MeshScalarT; - - // Input: - PHX::MDField wBF; - PHX::MDField wGradBF; - - PHX::MDField qFluct; //vector q' containing fluid fluctuations in primitive variables - PHX::MDField qFluctGrad; - PHX::MDField qFluctDot; - PHX::MDField force; - - PHX::MDField mu; - PHX::MDField lambda; - PHX::MDField kappa; - PHX::MDField tau11; - PHX::MDField tau12; - PHX::MDField tau13; - PHX::MDField tau22; - PHX::MDField tau23; - PHX::MDField tau33; - - double gamma_gas; //1.4 typically - double Rgas; //Non-dimensional gas constant Rgas = R*Tref/(cref*cref), where R = nondimensional gas constant = 287.0 typically - double Re; //Reynolds number - double Pr; //Prandtl number, 0.72 typically - - // Output: - PHX::MDField Residual; - - - std::size_t numNodes; - std::size_t numQPs; - std::size_t numDims; - std::size_t vecDim; - bool enableTransient; -}; -} - -#endif diff --git a/src/evaluators/pde/PHAL_ComprNSResid_Def.hpp b/src/evaluators/pde/PHAL_ComprNSResid_Def.hpp deleted file mode 100644 index 978014edad..0000000000 --- a/src/evaluators/pde/PHAL_ComprNSResid_Def.hpp +++ /dev/null @@ -1,202 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - - -#include "Teuchos_TestForException.hpp" -#include "Phalanx_DataLayout.hpp" - -#include "Intrepid2_FunctionSpaceTools.hpp" - -namespace PHAL { - - -//********************************************************************** -template -ComprNSResid:: -ComprNSResid(const Teuchos::ParameterList& p) : - wBF (p.get ("Weighted BF Name"), - p.get >("Node QP Scalar Data Layout") ), - wGradBF (p.get ("Weighted Gradient BF Name"), - p.get >("Node QP Gradient Data Layout") ), - qFluct (p.get ("QP Variable Name"), - p.get >("QP Vector Data Layout") ), - qFluctGrad (p.get ("Gradient QP Variable Name"), - p.get >("QP Tensor Data Layout") ), - qFluctDot (p.get ("QP Time Derivative Variable Name"), - p.get >("QP Vector Data Layout") ), - force (p.get ("Body Force Name"), - p.get >("QP Vector Data Layout") ), - mu (p.get ("Viscosity Mu QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - lambda (p.get ("Viscosity Lambda QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - kappa (p.get ("Viscosity Kappa QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau11 (p.get ("Stress Tau11 QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau12 (p.get ("Stress Tau12 QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau13 (p.get ("Stress Tau13 QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau22 (p.get ("Stress Tau22 QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau23 (p.get ("Stress Tau23 QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau33 (p.get ("Stress Tau33 QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - Residual (p.get ("Residual Name"), - p.get >("Node Vector Data Layout") ) -{ - - //TO DOs: - //3D - - this->addDependentField(qFluct.fieldTag()); - this->addDependentField(qFluctGrad.fieldTag()); - this->addDependentField(qFluctDot.fieldTag()); - this->addDependentField(force.fieldTag()); - this->addDependentField(mu.fieldTag()); - this->addDependentField(kappa.fieldTag()); - this->addDependentField(lambda.fieldTag()); - this->addDependentField(wBF.fieldTag()); - this->addDependentField(wGradBF.fieldTag()); - this->addDependentField(tau11.fieldTag()); - this->addDependentField(tau12.fieldTag()); - this->addDependentField(tau13.fieldTag()); - this->addDependentField(tau22.fieldTag()); - this->addDependentField(tau23.fieldTag()); - this->addDependentField(tau33.fieldTag()); - - this->addEvaluatedField(Residual); - - - this->setName("ComprNSResid" ); - - std::vector dims; - wGradBF.fieldTag().dataLayout().dimensions(dims); - numNodes = dims[1]; - numQPs = dims[2]; - numDims = dims[3]; - - - Teuchos::ParameterList* bf_list = - p.get("Parameter List"); - - qFluct.fieldTag().dataLayout().dimensions(dims); - vecDim = dims[2]; - - gamma_gas = bf_list->get("Gamma", 1.4); - Rgas = bf_list->get("Gas constant R", 0.714285733); - Pr = bf_list->get("Prandtl number Pr", 0.72); - Re = bf_list->get("Reynolds number Re", 1.0); - - - -std::cout << " vecDim = " << vecDim << std::endl; -std::cout << " numDims = " << numDims << std::endl; - - -if (vecDim != numDims+2) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, - std::endl << "Error in PHAL::ComprNS constructor: " << - "Invalid Parameter vecDim. vecDim should be numDims + 2 = " << numDims + 2 << "." << std::endl);} - - -} - -//********************************************************************** -template -void ComprNSResid:: -postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& fm) -{ - this->utils.setFieldData(qFluct,fm); - this->utils.setFieldData(qFluctGrad,fm); - this->utils.setFieldData(qFluctDot,fm); - this->utils.setFieldData(force,fm); - this->utils.setFieldData(mu,fm); - this->utils.setFieldData(kappa,fm); - this->utils.setFieldData(lambda,fm); - this->utils.setFieldData(wBF,fm); - this->utils.setFieldData(wGradBF,fm); - this->utils.setFieldData(tau11,fm); - this->utils.setFieldData(tau12,fm); - this->utils.setFieldData(tau13,fm); - this->utils.setFieldData(tau22,fm); - this->utils.setFieldData(tau23,fm); - this->utils.setFieldData(tau33,fm); - - this->utils.setFieldData(Residual,fm); -} - -//********************************************************************** -template -void ComprNSResid:: -evaluateFields(typename Traits::EvalData workset) -{ - typedef Intrepid2::FunctionSpaceTools FST; - - if (numDims == 2) { //2D case; order of variables is rho, u, v, T - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t node=0; node < numNodes; ++node) { - for (std::size_t i=0; i -class ComprNSViscosity : public PHX::EvaluatorWithBaseImpl, - public PHX::EvaluatorDerived { - -public: - - typedef typename EvalT::ScalarT ScalarT; - - ComprNSViscosity(const Teuchos::ParameterList& p); - - void postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& vm); - - void evaluateFields(typename Traits::EvalData d); - - -private: - - typedef typename EvalT::MeshScalarT MeshScalarT; - - // Input: - PHX::MDField coordVec; - PHX::MDField qFluct; //vector q' containing fluid fluctuations in primitive variables - PHX::MDField qFluctGrad; - //reference values for viscosities - double muref; - double kapparef; - double Tref; //reference temperature -- needed for Sutherland's viscosity law - double Pr; //Prandtl number - double Cp; //specific heat at constant pressure - - // Output: - PHX::MDField mu; - PHX::MDField kappa; - PHX::MDField lambda; - PHX::MDField tau11; - PHX::MDField tau12; - PHX::MDField tau13; - PHX::MDField tau22; - PHX::MDField tau23; - PHX::MDField tau33; - - //Force types - enum VISCTYPE {CONSTANT, SUTHERLAND}; - VISCTYPE visc_type; - - std::size_t numQPs; - std::size_t numDims; - std::size_t vecDim; - -}; -} - -#endif diff --git a/src/evaluators/pde/PHAL_ComprNSViscosity_Def.hpp b/src/evaluators/pde/PHAL_ComprNSViscosity_Def.hpp deleted file mode 100644 index 0b5397d1b6..0000000000 --- a/src/evaluators/pde/PHAL_ComprNSViscosity_Def.hpp +++ /dev/null @@ -1,178 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Teuchos_TestForException.hpp" -#include "Phalanx_DataLayout.hpp" -#include "Sacado.hpp" - -#include "Intrepid2_FunctionSpaceTools.hpp" - -namespace PHAL { -//********************************************************************** - -template -ComprNSViscosity:: -ComprNSViscosity(const Teuchos::ParameterList& p) : - qFluct (p.get ("QP Variable Name"), - p.get >("QP Vector Data Layout") ), - qFluctGrad (p.get ("Gradient QP Variable Name"), - p.get >("QP Tensor Data Layout") ), - mu (p.get ("Viscosity Mu QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - kappa (p.get ("Viscosity Kappa QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - lambda (p.get ("Viscosity Lambda QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau11 (p.get ("Stress Tau11 QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau12 (p.get ("Stress Tau12 QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau13 (p.get ("Stress Tau13 QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau22 (p.get ("Stress Tau22 QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau23 (p.get ("Stress Tau23 QP Variable Name"), - p.get >("QP Scalar Data Layout") ), - tau33 (p.get ("Stress Tau33 QP Variable Name"), - p.get >("QP Scalar Data Layout") ) -{ - Teuchos::ParameterList* visc_list = - p.get("Parameter List"); - - std::string viscType = visc_list->get("Type", "Constant"); - - if (viscType == "Constant"){ - std::cout << "Constant viscosity!" << std::endl; - visc_type = CONSTANT; - } - else if (viscType == "Sutherland") { - std::cout << "Sutherland viscosity!" << std::endl; - visc_type = SUTHERLAND; - } - - muref = visc_list->get("Mu_ref", 1.0); - kapparef = visc_list->get("Kappa_ref", 1.0); - Tref = visc_list->get("T_ref", 1.0); - Pr = visc_list->get("Prandtl number Pr", 0.72); - Cp = visc_list->get("Specific heat Cp", 1.0); - - coordVec = decltype(coordVec)( - p.get("Coordinate Vector Name"), - p.get >("QP Gradient Data Layout") ); - - this->addDependentField(qFluct.fieldTag()); - this->addDependentField(qFluctGrad.fieldTag()); - this->addDependentField(coordVec.fieldTag()); - this->addEvaluatedField(mu); - this->addEvaluatedField(kappa); - this->addEvaluatedField(lambda); - this->addEvaluatedField(tau11); - this->addEvaluatedField(tau12); - this->addEvaluatedField(tau13); - this->addEvaluatedField(tau22); - this->addEvaluatedField(tau23); - this->addEvaluatedField(tau33); - - std::vector dims; - qFluctGrad.fieldTag().dataLayout().dimensions(dims); - numQPs = dims[2]; - numDims = dims[3]; - - qFluct.fieldTag().dataLayout().dimensions(dims); - vecDim = dims[2]; - - std::cout << "vecdim in viscosity evaluator: " << vecDim << std::endl; - std::cout << "numDims in viscosity evaluator: " << numDims << std::endl; - std::cout << "numQPs in viscosity evaluator: " << numQPs << std::endl; - std::cout << "Mu_ref: " << muref << std::endl; - std::cout << "Kappa_ref: " << kapparef << std::endl; - - this->setName("ComprNSViscosity" ); -} - -//********************************************************************** -template -void ComprNSViscosity:: -postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& fm) -{ - this->utils.setFieldData(qFluct,fm); - this->utils.setFieldData(qFluctGrad,fm); - this->utils.setFieldData(mu,fm); - this->utils.setFieldData(kappa,fm); - this->utils.setFieldData(lambda,fm); - this->utils.setFieldData(coordVec,fm); - this->utils.setFieldData(tau11,fm); - this->utils.setFieldData(tau12,fm); - this->utils.setFieldData(tau13,fm); - this->utils.setFieldData(tau22,fm); - this->utils.setFieldData(tau23,fm); - this->utils.setFieldData(tau33,fm); -} - -//********************************************************************** -template -void ComprNSViscosity:: -evaluateFields(typename Traits::EvalData workset) -{ - //Visocisity coefficients - if (visc_type == CONSTANT){ - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t qp=0; qp < numQPs; ++qp) { - mu(cell,qp) = 1.0; - kappa(cell,qp) = mu(cell,qp)*Cp/Pr/kapparef; - mu(cell,qp) = 1.0/muref; //non-dimensionalize mu - lambda(cell,qp) = -2.0/3.0*mu(cell,qp); //Stokes' hypothesis - } - } - } - else if (visc_type == SUTHERLAND){ - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t qp=0; qp < numQPs; ++qp) { - ScalarT T = qFluct(cell,qp,vecDim-1)*Tref; //temperature (dimensional) - mu(cell,qp) = (1.458e-6)*sqrt(T*T*T)/(T + 110.4); //mu = (1.458e-6)*T^(1/5)/(T + 110.4) - kappa(cell,qp) = mu(cell,qp)*Cp/Pr/kapparef; - mu(cell,qp) = mu(cell,qp)/muref; //non-dimensionalize mu - lambda(cell,qp) = -2.0/3.0*mu(cell,qp); //Stokes' hypothesis - } - } - } - //Viscous stresses - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t qp=0; qp < numQPs; ++qp) { - tau11(cell,qp) = mu(cell,qp)*2.0*qFluctGrad(cell,qp,1,0) + lambda(cell,qp)*(qFluctGrad(cell,qp,1,0) + qFluctGrad(cell,qp,2,1)); //mu*2*du/dx + lambda*div(u) - tau12(cell,qp) = mu(cell,qp)*(qFluctGrad(cell,qp,1,1) + qFluctGrad(cell,qp,2,0)); //mu*(du/dy + dv/dx) - tau13(cell,qp) = 0.0; - tau22(cell,qp) = mu(cell,qp)*2.0*qFluctGrad(cell,qp,2,1) + lambda(cell,qp)*(qFluctGrad(cell,qp,1,0) + qFluctGrad(cell,qp,2,1)); //mu*2*dv/dy + lambda*div(u) - tau23(cell,qp) = 0.0; - tau33(cell,qp) = 0.0; - } - } - if (numDims == 3) {//3D case - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t qp=0; qp < numQPs; ++qp) { - tau11(cell,qp) += lambda(cell,qp)*qFluctGrad(cell,qp,3,2); //+lambda*dw/dz - tau13(cell,qp) += mu(cell,qp)*(qFluctGrad(cell,qp,1,2) + qFluctGrad(cell,qp,3,0)); //mu*(du/dz + dw/dx) - tau22(cell,qp) += lambda(cell,qp)*qFluctGrad(cell,qp,3,2); //+lambda*dw/dz - tau23(cell,qp) += mu(cell,qp)*(qFluctGrad(cell,qp,2,3) + qFluctGrad(cell,qp,3,1)); //mu*(dv/dz + dw/dy) - TEUCHOS_TEST_FOR_EXCEPTION( - true, std::logic_error, - "This next line has qFluct in it with the wrong indexing: there" - " should be 3, not 4. Inspection does not reveal what should be" - " fixed. I suspect qFluct should be qFluctGrad, but I can't be" - " sure. I suspect there is no test coverage of this codepath, so" - " for now I'll do the safe thing and throw an exception. I also" - " have to inactivate the code, as it won't compile with Kokkos."); -#if 0 - tau33(cell,qp) += 2.0*mu(cell,qp)*qFluctGrad(cell,qp,3,2) + lambda(cell,qp)*(qFluctGrad(cell,qp,1,0) + qFluctGrad(cell,qp,2,1) + qFluct(cell,qp,3,2)); //mu*2*dw/dz + lambda*div(u) -#endif - } - } - } -} - -} - diff --git a/src/evaluators/pde/PHAL_LinComprNSBodyForce.cpp b/src/evaluators/pde/PHAL_LinComprNSBodyForce.cpp deleted file mode 100644 index e9706186e2..0000000000 --- a/src/evaluators/pde/PHAL_LinComprNSBodyForce.cpp +++ /dev/null @@ -1,13 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "PHAL_AlbanyTraits.hpp" - -#include "PHAL_LinComprNSBodyForce.hpp" -#include "PHAL_LinComprNSBodyForce_Def.hpp" - -PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::LinComprNSBodyForce) - diff --git a/src/evaluators/pde/PHAL_LinComprNSBodyForce.hpp b/src/evaluators/pde/PHAL_LinComprNSBodyForce.hpp deleted file mode 100644 index 53d83ac9ad..0000000000 --- a/src/evaluators/pde/PHAL_LinComprNSBodyForce.hpp +++ /dev/null @@ -1,60 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef PHAL_LINCOMPRNSBODYFORCE_HPP -#define PHAL_LINCOMPRNSBODYFORCE_HPP - -#include "Phalanx_config.hpp" -#include "Phalanx_Evaluator_WithBaseImpl.hpp" -#include "Phalanx_Evaluator_Derived.hpp" -#include "Phalanx_MDField.hpp" - -namespace PHAL { -/** \brief Finite Element Interpolation Evaluator - - This evaluator interpolates nodal DOF values to quad points. - -*/ - -template -class LinComprNSBodyForce : public PHX::EvaluatorWithBaseImpl, - public PHX::EvaluatorDerived { - -public: - - typedef typename EvalT::ScalarT ScalarT; - - LinComprNSBodyForce(const Teuchos::ParameterList& p); - - void postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& vm); - - void evaluateFields(typename Traits::EvalData d); - - -private: - - typedef typename EvalT::MeshScalarT MeshScalarT; - - // Input: - PHX::MDField coordVec; - Teuchos::Array gravity; - - // Output: - PHX::MDField force; - - //Force types - enum BFTYPE {NONE, STEADYEUL, UNSTEADYEULMMS, DRIVENPULSE}; - BFTYPE bf_type; - - std::size_t numQPs; - std::size_t numDims; - std::size_t vecDim; - -}; -} - -#endif diff --git a/src/evaluators/pde/PHAL_LinComprNSBodyForce_Def.hpp b/src/evaluators/pde/PHAL_LinComprNSBodyForce_Def.hpp deleted file mode 100644 index 13f0033ee6..0000000000 --- a/src/evaluators/pde/PHAL_LinComprNSBodyForce_Def.hpp +++ /dev/null @@ -1,143 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Teuchos_TestForException.hpp" -#include "Phalanx_DataLayout.hpp" -#include "Sacado.hpp" - -#include "Intrepid2_FunctionSpaceTools.hpp" - -namespace PHAL { -const double pi = 3.1415926535897932385; -//********************************************************************** - -template -LinComprNSBodyForce:: -LinComprNSBodyForce(const Teuchos::ParameterList& p) : - force(p.get("Body Force Name"), - p.get >("QP Vector Data Layout") ) -{ - std::cout << "Lin Compr NS body force constructor!" << std::endl; - Teuchos::ParameterList* bf_list = - p.get("Parameter List"); - - std::string type = bf_list->get("Type", "None"); - if (type == "None") bf_type = NONE; - else if (type == "Steady Euler") bf_type = STEADYEUL; - else if (type == "Unsteady Euler MMS") bf_type = UNSTEADYEULMMS; - else if (type == "Driven Pulse") bf_type = DRIVENPULSE; - - if (bf_type != NONE) { - coordVec = decltype(coordVec)( - p.get("Coordinate Vector Name"), - p.get >("QP Gradient Data Layout") ); - this->addDependentField(coordVec); - } - - this->addEvaluatedField(force); - - Teuchos::RCP gradient_dl = - p.get< Teuchos::RCP >("QP Gradient Data Layout"); - std::vector dims; - gradient_dl->dimensions(dims); - numQPs = dims[1]; - numDims = dims[2]; - Teuchos::RCP vector_dl = - p.get< Teuchos::RCP >("QP Vector Data Layout"); - vector_dl->dimensions(dims); - vecDim = dims[2]; - -std::cout << " in Lin Compr NS Stokes source! " << std::endl; -std::cout << " vecDim = " << vecDim << std::endl; -std::cout << " numDims = " << numDims << std::endl; -std::cout << " numQPs = " << numQPs << std::endl; - - - this->setName("LinComprNSBodyForce" ); -} - -//********************************************************************** -template -void LinComprNSBodyForce:: -postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& fm) -{ - if (bf_type == STEADYEUL || bf_type == UNSTEADYEULMMS || bf_type == DRIVENPULSE) { - this->utils.setFieldData(coordVec,fm); - } - - this->utils.setFieldData(force,fm); -} - -//********************************************************************** -template -void LinComprNSBodyForce:: -evaluateFields(typename Traits::EvalData workset) -{ - if (bf_type == NONE) { - for (std::size_t cell=0; cell < workset.numCells; ++cell) - for (std::size_t qp=0; qp < numQPs; ++qp) - for (std::size_t i=0; i < vecDim; ++i) - force(cell,qp,i) = 0.0; - } - else if (bf_type == STEADYEUL) { - const double ubar = 1.0; - const double vbar = 1.0; - const double zetabar = 1.0; - const double pbar = 0.0; - const double gamma_gas = 1.4; - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t qp=0; qp < numQPs; ++qp) { - MeshScalarT x = coordVec(cell,qp,0); - MeshScalarT y = coordVec(cell,qp,1); - force(cell,qp,0) = -1.0*(ubar*(y - x*sin(x)) + vbar*x + zetabar*2.0*x*(0.5-y)); - force(cell,qp,1) = -1.0*(ubar*cos(x)*y + vbar*sin(x) - zetabar*x*x); - force(cell,qp,2) = -1.0*(gamma_gas*pbar*(y - x*sin(x) + sin(x)) + ubar*2.0*x*(0.5-y) - vbar*x*x); - } - } - } - else if (bf_type == UNSTEADYEULMMS) { - const double ubar = 0.0; - const double vbar = 0.0; - const double zetabar = 1.0; - const double pbar = 0.7142857; - const double a = 1.0; - const double gamma_gas = 1.4; - const RealType time = workset.current_time; - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t qp=0; qp < numQPs; ++qp) { - MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0); - MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1); - force(cell,qp,0) = -1.0*exp(-a*time)*(-a*sin(x2pi)*cos(y2pi) + ubar*2.0*pi*cos(x2pi)*cos(y2pi) - -vbar*2.0*pi*sin(x2pi)*sin(y2pi) + 2.0*pi*zetabar*cos(x2pi)*sin(y2pi)); - force(cell,qp,1) = -1.0*exp(-a*time)*(-a*cos(x2pi)*sin(y2pi) - 2.0*pi*ubar*sin(x2pi)*sin(y2pi) - + vbar*2.0*pi*cos(x2pi)*cos(y2pi) + 2.0*pi*zetabar*sin(x2pi)*cos(y2pi)); - force(cell,qp,2) = -1.0*exp(-a*time)*(-a*sin(x2pi)*sin(y2pi) + gamma_gas*pbar*4.0*pi*cos(x2pi)*cos(y2pi) + - ubar*2.0*pi*cos(x2pi)*sin(y2pi) + vbar*2.0*pi*sin(x2pi)*cos(y2pi)); - } - } - } - else if (bf_type == DRIVENPULSE) { - const RealType time = workset.current_time; - const double tref = 1.0/347.9693; - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t qp=0; qp < numQPs; ++qp) { - MeshScalarT x = coordVec(cell,qp,0); - MeshScalarT y = coordVec(cell,qp,1); - force(cell,qp,0) = 0.0; - if ((x >= 0.9) && (x <= 1.0) && (y >= 0.9) && (y <= 1.0)) - force(cell,qp,1) = (1.0e-4)*cos(2.0*pi*1000*time*tref); - else - force(cell,qp,1) = 0.0; - force(cell,qp,2) = 0.0; - } - } - - } -} - -} - diff --git a/src/evaluators/pde/PHAL_LinComprNSResid.cpp b/src/evaluators/pde/PHAL_LinComprNSResid.cpp deleted file mode 100644 index 5a60b84a28..0000000000 --- a/src/evaluators/pde/PHAL_LinComprNSResid.cpp +++ /dev/null @@ -1,13 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "PHAL_AlbanyTraits.hpp" - -#include "PHAL_LinComprNSResid.hpp" -#include "PHAL_LinComprNSResid_Def.hpp" - -PHAL_INSTANTIATE_TEMPLATE_CLASS(PHAL::LinComprNSResid) - diff --git a/src/evaluators/pde/PHAL_LinComprNSResid.hpp b/src/evaluators/pde/PHAL_LinComprNSResid.hpp deleted file mode 100644 index 038b2b58b5..0000000000 --- a/src/evaluators/pde/PHAL_LinComprNSResid.hpp +++ /dev/null @@ -1,72 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef PHAL_LINCOMPRNSRESID_HPP -#define PHAL_LINCOMPRNSRESID_HPP - -#include "Phalanx_config.hpp" -#include "Phalanx_Evaluator_WithBaseImpl.hpp" -#include "Phalanx_Evaluator_Derived.hpp" -#include "Phalanx_MDField.hpp" - -namespace PHAL { -/** \brief Finite Element Interpolation Evaluator - - This evaluator interpolates nodal DOF values to quad points. - -*/ - -template -class LinComprNSResid : public PHX::EvaluatorWithBaseImpl, - public PHX::EvaluatorDerived { - -public: - - LinComprNSResid(const Teuchos::ParameterList& p); - - void postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& vm); - - void evaluateFields(typename Traits::EvalData d); - -private: - - typedef typename EvalT::ScalarT ScalarT; - typedef typename EvalT::MeshScalarT MeshScalarT; - - // Input: - PHX::MDField wBF; - PHX::MDField wGradBF; - - PHX::MDField qFluct; //vector q' containing fluid fluctuations in primitive variables - PHX::MDField qFluctGrad; - PHX::MDField qFluctDot; - PHX::MDField force; - - Teuchos::Array baseFlowData; - double gamma_gas; //1.4 typically - double Rgas; //Non-dimensional gas constant Rgas = R*Tref/(cref*cref), where R = nondimensional gas constant = 287.0 typically - double Re; //Reynolds number - double Pr; //Prandtl number, 0.72 typically - double mu; double lambda; //viscosity coefficients - double kappa; //thermal diffusivity - bool IBP_convect_terms; //boolean specifying whether you want to integrate by parts the convective terms in the weak form - - // Output: - PHX::MDField Residual; - - - std::size_t numNodes; - std::size_t numQPs; - std::size_t numDims; - std::size_t vecDim; - bool enableTransient; - enum EQNTYPE {EULER, NS}; - EQNTYPE eqn_type; -}; -} - -#endif diff --git a/src/evaluators/pde/PHAL_LinComprNSResid_Def.hpp b/src/evaluators/pde/PHAL_LinComprNSResid_Def.hpp deleted file mode 100644 index b009458162..0000000000 --- a/src/evaluators/pde/PHAL_LinComprNSResid_Def.hpp +++ /dev/null @@ -1,422 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - - -#include "Teuchos_TestForException.hpp" -#include "Phalanx_DataLayout.hpp" - -#include "Intrepid2_FunctionSpaceTools.hpp" - -namespace PHAL { - - -//********************************************************************** -template -LinComprNSResid:: -LinComprNSResid(const Teuchos::ParameterList& p) : - wBF (p.get ("Weighted BF Name"), - p.get >("Node QP Scalar Data Layout") ), - wGradBF (p.get ("Weighted Gradient BF Name"), - p.get >("Node QP Gradient Data Layout") ), - qFluct (p.get ("QP Variable Name"), - p.get >("QP Vector Data Layout") ), - qFluctGrad (p.get ("Gradient QP Variable Name"), - p.get >("QP Tensor Data Layout") ), - qFluctDot (p.get ("QP Time Derivative Variable Name"), - p.get >("QP Vector Data Layout") ), - force (p.get ("Body Force Name"), - p.get >("QP Vector Data Layout") ), - Residual (p.get ("Residual Name"), - p.get >("Node Vector Data Layout") ) -{ - - - if (p.isType("Disable Transient")) - enableTransient = !p.get("Disable Transient"); - else enableTransient = true; - - this->addDependentField(qFluct.fieldTag()); - this->addDependentField(qFluctGrad.fieldTag()); - if(enableTransient) - this->addDependentField(qFluctDot.fieldTag()); - this->addDependentField(force.fieldTag()); - this->addDependentField(wBF.fieldTag()); - this->addDependentField(wGradBF.fieldTag()); - - this->addEvaluatedField(Residual); - - - this->setName("LinComprNSResid" ); - - std::vector dims; - wGradBF.fieldTag().dataLayout().dimensions(dims); - numNodes = dims[1]; - numQPs = dims[2]; - numDims = dims[3]; - - - Teuchos::ParameterList* bf_list = - p.get("Parameter List"); - std::string eqnType = bf_list->get("Type", "Euler"); - - if (eqnType == "Euler") { - std::cout << "setting euler equations!" << std::endl; - eqn_type = EULER; - } - else if (eqnType == "Navier-Stokes") { - std::cout << "setting n-s equations!" << std::endl; - eqn_type = NS; - } - - - qFluct.fieldTag().dataLayout().dimensions(dims); - vecDim = dims[2]; - - Teuchos::Array defaultBaseFlowData(numDims+2); - baseFlowData = bf_list->get("Base Flow Data", defaultBaseFlowData); - //for EULER, baseFlowData = (ubar, vbar, wbar, zetabar, pbar) - //for NS, baseFlowData = (ubar, vbar, wbar, Tbar, rhobar) - - gamma_gas = bf_list->get("Gamma", 1.4); - Rgas = bf_list->get("Gas constant R", 0.714285733); - Pr = bf_list->get("Prandtl number Pr", 0.72); - Re = bf_list->get("Reynolds number Re", 1.0); - mu = bf_list->get("Viscocity mu", 0.0); - lambda = -2.0/3.0*mu; //Stokes' hypothesis - kappa = bf_list->get("Diffusivity kappa", 0.0); - IBP_convect_terms = bf_list->get("IBP Convective Terms", false); - - if (IBP_convect_terms == true) - std::cout << "Integrating convective terms by parts in weak form." << std::endl; - - -std::cout << " vecDim = " << vecDim << std::endl; -std::cout << " numDims = " << numDims << std::endl; - -if (baseFlowData.size()!=static_cast(numDims+2)) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, - std::endl << "Error in PHAL::LinComprNS constructor: " << - "baseFlow data should have length numDims + 2 = " << numDims+2 << "." << std::endl);} - - -if ((eqn_type == EULER) && (vecDim != numDims+1)) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, - std::endl << "Error in PHAL::LinComprNS constructor: " << - "Invalid Parameter vecDim. vecDim should be numDims + 1 = " << numDims + 1 << " for Euler equations." << std::endl);} - -if ((eqn_type == NS) && (vecDim != numDims+2)) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, - std::endl << "Error in PHAL::LinComprNS constructor: " << - "Invalid Parameter vecDim. vecDim should be numDims + 2 = " << numDims + 2 << " for Navier-Stokes equations." << std::endl);} - -} - -//********************************************************************** -template -void LinComprNSResid:: -postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& fm) -{ - this->utils.setFieldData(qFluct,fm); - this->utils.setFieldData(qFluctGrad,fm); - if(enableTransient) - this->utils.setFieldData(qFluctDot,fm); - this->utils.setFieldData(force,fm); - this->utils.setFieldData(wBF,fm); - this->utils.setFieldData(wGradBF,fm); - - this->utils.setFieldData(Residual,fm); -} - -//********************************************************************** -template -void LinComprNSResid:: -evaluateFields(typename Traits::EvalData workset) -{ - if (eqn_type == EULER) { //Euler equations - if (numDims == 1) { //1D case - double ubar = baseFlowData[0]; - double zetabar = baseFlowData[1]; - double pbar = baseFlowData[2]; - if (IBP_convect_terms == false) {//variational formulation in which the convective terms are not integrated by parts - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t node=0; node < numNodes; ++node) { - for (std::size_t i=0; i -class ConcentrationResid : public PHX::EvaluatorWithBaseImpl, - public PHX::EvaluatorDerived { - -public: - - ConcentrationResid(const Teuchos::ParameterList& p, - const Teuchos::RCP& dl); - - void postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& vm); - - void evaluateFields(typename Traits::EvalData d); - -private: - - typedef typename EvalT::ScalarT ScalarT; - typedef typename EvalT::MeshScalarT MeshScalarT; - - // Input: - PHX::MDField wBF; - PHX::MDField wGradBF; - PHX::MDField PotentialGrad; - PHX::MDField Concentration; - PHX::MDField Concentration_dot; - PHX::MDField ConcentrationGrad; - - // Output: - PHX::MDField ConcentrationResidual; - - unsigned int numNodes, numQPs, numDims, numSpecies; - std::vector D,beta; // Placeholder for charges - - bool enableTransient; -}; -} - -#endif diff --git a/src/evaluators/pde/PNP_ConcentrationResid_Def.hpp b/src/evaluators/pde/PNP_ConcentrationResid_Def.hpp deleted file mode 100644 index 3ce216cedb..0000000000 --- a/src/evaluators/pde/PNP_ConcentrationResid_Def.hpp +++ /dev/null @@ -1,104 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Teuchos_TestForException.hpp" -#include "Phalanx_DataLayout.hpp" - -#include "Intrepid2_FunctionSpaceTools.hpp" - - -//********************************************************************** -template -PNP::ConcentrationResid:: -ConcentrationResid(const Teuchos::ParameterList& p, - const Teuchos::RCP& dl) : - wBF (p.get ("Weighted BF Name"), dl->node_qp_scalar), - wGradBF (p.get ("Weighted Gradient BF Name"), dl->node_qp_gradient), - PotentialGrad ("Potential Gradient", dl->qp_gradient), - Concentration ("Concentration", dl->qp_vector), - Concentration_dot ("Concentration_dot", dl->qp_vector), - ConcentrationGrad ("Concentration Gradient", dl->qp_vecgradient), - ConcentrationResidual ("Concentration Residual", dl->node_vector ) -{ - if (p.isType("Disable Transient")) - enableTransient = !p.get("Disable Transient"); - else enableTransient = true; - - this->addDependentField(wBF.fieldTag()); - this->addDependentField(wGradBF.fieldTag()); - this->addDependentField(Concentration.fieldTag()); - if (enableTransient) this->addDependentField(Concentration_dot.fieldTag()); - this->addDependentField(ConcentrationGrad.fieldTag()); - this->addDependentField(PotentialGrad.fieldTag()); - - this->addEvaluatedField(ConcentrationResidual); - - std::vector dims; - wGradBF.fieldTag().dataLayout().dimensions(dims); - numNodes = dims[1]; - numQPs = dims[2]; - numDims = dims[3]; - ConcentrationGrad.fieldTag().dataLayout().dimensions(dims); - numSpecies = dims[2]; - - // Placeholder for properties - beta.resize(numSpecies); - beta[0] = 1.0; - beta[1] = -1.0; - D.resize(numSpecies); - D[0] = 1.0; - D[1] = 2.0; - - this->setName("ConcentrationResid" ); -} - -//********************************************************************** -template -void PNP::ConcentrationResid:: -postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& fm) -{ - this->utils.setFieldData(wBF,fm); - this->utils.setFieldData(wGradBF,fm); - this->utils.setFieldData(Concentration,fm); - if (enableTransient) this->utils.setFieldData(Concentration_dot,fm); - this->utils.setFieldData(ConcentrationGrad,fm); - this->utils.setFieldData(PotentialGrad,fm); - - this->utils.setFieldData(ConcentrationResidual,fm); -} - -//********************************************************************** -template -void PNP::ConcentrationResid:: -evaluateFields(typename Traits::EvalData workset) -{ - typedef Intrepid2::FunctionSpaceTools FST; - - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t node=0; node < numNodes; ++node) { - for (std::size_t j=0; j < numSpecies; ++j) { - ConcentrationResidual(cell,node,j) = 0.0; - } } } - - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t node=0; node < numNodes; ++node) { - for (std::size_t qp=0; qp < numQPs; ++qp) { - for (std::size_t j=0; j < numSpecies; ++j) { - for (std::size_t dim=0; dim < numDims; ++dim) { - ConcentrationResidual(cell,node,j) += - D[j]*(ConcentrationGrad(cell,qp,j,dim) - + beta[j]*Concentration(cell,qp,j)*PotentialGrad(cell,qp,dim)) - *wGradBF(cell,node,qp,dim); - } - } - } - } - } - -} -//********************************************************************** - diff --git a/src/evaluators/pde/PNP_PotentialResid.cpp b/src/evaluators/pde/PNP_PotentialResid.cpp deleted file mode 100644 index ed87066dd0..0000000000 --- a/src/evaluators/pde/PNP_PotentialResid.cpp +++ /dev/null @@ -1,13 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "PHAL_AlbanyTraits.hpp" - -#include "PNP_PotentialResid.hpp" -#include "PNP_PotentialResid_Def.hpp" - -PHAL_INSTANTIATE_TEMPLATE_CLASS(PNP::PotentialResid) - diff --git a/src/evaluators/pde/PNP_PotentialResid.hpp b/src/evaluators/pde/PNP_PotentialResid.hpp deleted file mode 100644 index 170e2e1f48..0000000000 --- a/src/evaluators/pde/PNP_PotentialResid.hpp +++ /dev/null @@ -1,58 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef PNP_POISSONRESID_HPP -#define PNP_POISSONRESID_HPP - -#include "Phalanx_config.hpp" -#include "Phalanx_Evaluator_WithBaseImpl.hpp" -#include "Phalanx_Evaluator_Derived.hpp" -#include "Phalanx_MDField.hpp" - -#include "Albany_Layouts.hpp" - -namespace PNP { -/** \brief Finite Element Interpolation Evaluator - - This evaluator interpolates nodal DOF values to quad points. - -*/ - -template -class PotentialResid : public PHX::EvaluatorWithBaseImpl, - public PHX::EvaluatorDerived { - -public: - - PotentialResid(const Teuchos::ParameterList& p, - const Teuchos::RCP& dl); - - void postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& vm); - - void evaluateFields(typename Traits::EvalData d); - -private: - - typedef typename EvalT::ScalarT ScalarT; - typedef typename EvalT::MeshScalarT MeshScalarT; - - // Input: - PHX::MDField wBF; - PHX::MDField wGradBF; - PHX::MDField PotentialGrad; - PHX::MDField Concentration; - PHX::MDField Permittivity; - - // Output: - PHX::MDField PotentialResidual; - - unsigned int numNodes, numQPs, numSpecies; - std::vector q; // Placeholder for charges -}; -} - -#endif diff --git a/src/evaluators/pde/PNP_PotentialResid_Def.hpp b/src/evaluators/pde/PNP_PotentialResid_Def.hpp deleted file mode 100644 index aef3320a18..0000000000 --- a/src/evaluators/pde/PNP_PotentialResid_Def.hpp +++ /dev/null @@ -1,92 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Teuchos_TestForException.hpp" -#include "Phalanx_DataLayout.hpp" - -#include "Intrepid2_FunctionSpaceTools.hpp" -#include "PHAL_Utilities.hpp" - - -//********************************************************************** -template -PNP::PotentialResid:: -PotentialResid(const Teuchos::ParameterList& p, - const Teuchos::RCP& dl) : - wBF (p.get ("Weighted BF Name"), dl->node_qp_scalar), - wGradBF (p.get ("Weighted Gradient BF Name"), dl->node_qp_gradient), - PotentialGrad ("Potential Gradient", dl->qp_gradient), - Concentration ("Concentration", dl->qp_vector), - Permittivity (p.get ("Permittivity Name"), dl->qp_scalar), - PotentialResidual ("Potential Residual", dl->node_scalar ) -{ - this->addDependentField(wBF.fieldTag()); - this->addDependentField(wGradBF.fieldTag()); - this->addDependentField(Permittivity.fieldTag()); - this->addDependentField(Concentration.fieldTag()); - this->addDependentField(PotentialGrad.fieldTag()); - - this->addEvaluatedField(PotentialResidual); - - std::vector dims; - wBF.fieldTag().dataLayout().dimensions(dims); - numNodes = dims[1]; - numQPs = dims[2]; - Concentration.fieldTag().dataLayout().dimensions(dims); - numSpecies = dims[2]; - - // Placeholder for charges - q.resize(numSpecies); - q[0] = 5.0; - q[1] = -5.0; - - this->setName("PotentialResid" ); -} - -//********************************************************************** -template -void PNP::PotentialResid:: -postRegistrationSetup(typename Traits::SetupData d, - PHX::FieldManager& fm) -{ - this->utils.setFieldData(wBF,fm); - this->utils.setFieldData(wGradBF,fm); - this->utils.setFieldData(Permittivity,fm); - this->utils.setFieldData(Concentration,fm); - this->utils.setFieldData(PotentialGrad,fm); - - this->utils.setFieldData(PotentialResidual,fm); -} - -//********************************************************************** -template -void PNP::PotentialResid:: -evaluateFields(typename Traits::EvalData workset) -{ - typedef Intrepid2::FunctionSpaceTools FST; - - // Scale gradient into a flux - // can't reuse memory, dependent fields must be const - auto flux = PHAL::create_copy("tmp_flux", PotentialGrad.get_view()); - FST::scalarMultiplyDataData (flux, Permittivity.get_view(), PotentialGrad.get_view()); - FST::integrate(PotentialResidual.get_view(), flux, wGradBF.get_view(), false); // "false" overwrites - - - for (std::size_t cell=0; cell < workset.numCells; ++cell) { - for (std::size_t node=0; node < numNodes; ++node) { - for (std::size_t qp=0; qp < numQPs; ++qp) { - for (std::size_t j=0; j < numSpecies; ++j) { - PotentialResidual(cell,node) -= - q[j]*Concentration(cell,qp,j)*wBF(cell,node,qp); -//cout << "XXX " << cell << " " << node << " " << qp << " " << j << " " << q[j] << " " << Concentration(cell,qp,j) << endl; - } - } - } - } - -} -//********************************************************************** - diff --git a/src/problems/Albany_CahnHillProblem.cpp b/src/problems/Albany_CahnHillProblem.cpp deleted file mode 100644 index 7eb1bcf214..0000000000 --- a/src/problems/Albany_CahnHillProblem.cpp +++ /dev/null @@ -1,105 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Albany_CahnHillProblem.hpp" - -#include "Intrepid2_DefaultCubatureFactory.hpp" -#include "Shards_CellTopology.hpp" -#include "PHAL_FactoryTraits.hpp" -#include "Albany_Utils.hpp" -#include "Albany_BCUtils.hpp" - - -namespace Albany { - -CahnHillProblem:: -CahnHillProblem (const Teuchos::RCP& params_, - const Teuchos::RCP& paramLib_, - const int numDim_, - const Teuchos::RCP >& comm_) : - AbstractProblem(params_, paramLib_, 2), - numDim(numDim_), - haveNoise(false), - comm(comm_), - use_sdbcs_(false) -{} - -CahnHillProblem:: -~CahnHillProblem() -{ -} - -void -CahnHillProblem:: -buildProblem( - Teuchos::ArrayRCP > meshSpecs, - StateManager& stateMgr) -{ - /* Construct All Phalanx Evaluators */ - TEUCHOS_TEST_FOR_EXCEPTION(meshSpecs.size()!=1,std::logic_error,"Problem supports one Material Block"); - - fm.resize(1); - fm[0] = Teuchos::rcp(new PHX::FieldManager); - buildEvaluators(*fm[0], *meshSpecs[0], stateMgr, BUILD_RESID_FM, - Teuchos::null); - - if(meshSpecs[0]->nsNames.size() > 0) // Build a nodeset evaluator if nodesets are present - - constructDirichletEvaluators(meshSpecs[0]->nsNames); - -} - -Teuchos::Array > -CahnHillProblem:: -buildEvaluators( - PHX::FieldManager& fm0, - const MeshSpecsStruct& meshSpecs, - StateManager& stateMgr, - FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList) -{ - // Call constructEvaluators(*rfm[0], *meshSpecs[0], stateMgr); - // for each EvalT in PHAL::AlbanyTraits::BEvalTypes - ConstructEvaluatorsOp op( - *this, fm0, meshSpecs, stateMgr, fmchoice, responseList); - Sacado::mpl::for_each fe(op); - return *op.tags; -} - -// Dirichlet BCs -void -CahnHillProblem::constructDirichletEvaluators(const std::vector& nodeSetIDs) -{ - // Construct BC evaluators for all node sets and names - std::vector bcNames(neq); - bcNames[0] = "rho"; - BCUtils bcUtils; - dfm = bcUtils.constructBCEvaluators(nodeSetIDs, bcNames, - this->params, this->paramLib); - use_sdbcs_ = bcUtils.useSDBCs(); - offsets_ = bcUtils.getOffsets(); - nodeSetIDs_ = bcUtils.getNodeSetIDs(); -} - -Teuchos::RCP -CahnHillProblem::getValidProblemParameters() const -{ - Teuchos::RCP validPL = - this->getGenericProblemParams("ValidCahnHillProblemParams"); - - Teuchos::Array defaultPeriod; - - validPL->set("b", 0.0, "b value in equation 1.1"); - validPL->set("gamma", 0.0, "gamma value in equation 2.2"); - validPL->set("Langevin Noise SD", 0.0, "Standard deviation of the Langevin noise to apply"); - validPL->set >("Langevin Noise Time Period", defaultPeriod, - "Time period to apply Langevin noise"); - validPL->set("Lump Mass", true, "Lump mass matrix in time derivative term"); - - return validPL; -} - -} // namespace Albany diff --git a/src/problems/Albany_CahnHillProblem.hpp b/src/problems/Albany_CahnHillProblem.hpp deleted file mode 100644 index bb7857a840..0000000000 --- a/src/problems/Albany_CahnHillProblem.hpp +++ /dev/null @@ -1,318 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef ALBANY_CAHNHILLPROBLEM_HPP -#define ALBANY_CAHNHILLPROBLEM_HPP - -#include "Teuchos_RCP.hpp" -#include "Teuchos_ParameterList.hpp" - -#include "Albany_AbstractProblem.hpp" - -#include "PHAL_Workset.hpp" -#include "PHAL_Dimension.hpp" -#include "Albany_ProblemUtils.hpp" - -namespace Albany { - - /*! - * \brief Abstract interface for representing a 1-D finite element - * problem. - */ - class CahnHillProblem : public AbstractProblem { - public: - - //! Default constructor - CahnHillProblem (const Teuchos::RCP& params, - const Teuchos::RCP& paramLib, - const int numDim_, - const Teuchos::RCP >& comm); - - //! Destructor - ~CahnHillProblem(); - - //! Return number of spatial dimensions - virtual int spatialDimension() const { return numDim; } - - //! Get boolean telling code if SDBCs are utilized - virtual bool useSDBCs() const {return use_sdbcs_; } - - //! Build the PDE instantiations, boundary conditions, and initial solution - virtual void buildProblem( - Teuchos::ArrayRCP > meshSpecs, - StateManager& stateMgr); - - // Build evaluators - virtual Teuchos::Array< Teuchos::RCP > - buildEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList); - - //! Each problem must generate it's list of valid parameters - Teuchos::RCP getValidProblemParameters() const; - - CahnHillProblem(const CahnHillProblem&) = delete; - CahnHillProblem& operator=(const CahnHillProblem&) = delete; - - //! Main problem setup routine. Not directly called, but indirectly by following functions - template - Teuchos::RCP - constructEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList); - - void constructDirichletEvaluators(const std::vector& nodeSetIDs); - - protected: - - int numDim; - - bool haveNoise; // Langevin noise present - - Teuchos::RCP > comm; - - Teuchos::RCP dl; - - /// Boolean marking whether SDBCs are used - bool use_sdbcs_; - }; - -} - -#include "Intrepid2_DefaultCubatureFactory.hpp" -#include "Shards_CellTopology.hpp" -#include "Albany_Utils.hpp" -#include "Albany_ProblemUtils.hpp" -#include "Albany_EvaluatorUtils.hpp" -#include "Albany_ResponseUtilities.hpp" - -#include "PHAL_CahnHillChemTerm.hpp" -#include "PHAL_LangevinNoiseTerm.hpp" -#include "PHAL_CahnHillRhoResid.hpp" -#include "PHAL_CahnHillWResid.hpp" - - -template -Teuchos::RCP -Albany::CahnHillProblem::constructEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fieldManagerChoice, - const Teuchos::RCP& responseList) -{ - using Teuchos::RCP; - using Teuchos::rcp; - using Teuchos::ParameterList; - using PHX::DataLayout; - using PHX::MDALayout; - using std::vector; - using std::string; - using PHAL::AlbanyTraits; - - // Problem is transient - TEUCHOS_TEST_FOR_EXCEPTION( - number_of_time_deriv != 1, - std::logic_error, - "Albany_CahnHillProblem must be defined as a transient calculation."); - - const CellTopologyData * const elem_top = &meshSpecs.ctd; - - RCP > - intrepidBasis = Albany::getIntrepid2Basis(*elem_top); - RCP cellType = rcp(new shards::CellTopology (elem_top)); - - - const int numNodes = intrepidBasis->getCardinality(); - const int worksetSize = meshSpecs.worksetSize; - - Intrepid2::DefaultCubatureFactory cubFactory; - RCP > cellCubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree); - - const int numQPtsCell = cellCubature->getNumPoints(); - const int numVertices = cellType->getNodeCount(); - - - *out << "Field Dimensions: Workset=" << worksetSize - << ", Vertices= " << numVertices - << ", Nodes= " << numNodes - << ", QuadPts= " << numQPtsCell - << ", Dim= " << numDim << std::endl; - - dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPtsCell,numDim)); - Albany::EvaluatorUtils evalUtils(dl); - - // Temporary variable used numerous times below - Teuchos::RCP > ev; - - Teuchos::ArrayRCP dof_names(neq); - dof_names[0] = "Rho"; // The concentration difference variable 0 \leq \rho \leq 1 - dof_names[1] = "W"; // The chemical potential difference variable - Teuchos::ArrayRCP dof_names_dot(neq); - dof_names_dot[0] = "rhoDot"; - dof_names_dot[1] = "wDot"; // not currently used - Teuchos::ArrayRCP resid_names(neq); - resid_names[0] = "Rho Residual"; - resid_names[1] = "W Residual"; - - fm0.template registerEvaluator - (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot)); - - fm0.template registerEvaluator - (evalUtils.constructScatterResidualEvaluator(false, resid_names)); - - fm0.template registerEvaluator - (evalUtils.constructGatherCoordinateVectorEvaluator()); - - fm0.template registerEvaluator - (evalUtils.constructMapToPhysicalFrameEvaluator( cellType, cellCubature)); - - fm0.template registerEvaluator - (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cellCubature)); - - for (unsigned int i=0; i - (evalUtils.constructDOFInterpolationEvaluator(dof_names[i],i)); - - fm0.template registerEvaluator - (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[i],i)); - - fm0.template registerEvaluator - (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[i],i)); - } - - - { // Form the Chemical Energy term in Eq. 2.2 - - RCP p = rcp(new ParameterList("Chem Energy Term")); - - p->set >("Parameter Library", paramLib); - - // b value in Equation 1.1 - p->set("b Value", params->get("b")); - - //Input - p->set("Rho QP Variable Name", "Rho"); - p->set("W QP Variable Name", "W"); - - p->set< RCP >("QP Scalar Data Layout", dl->qp_scalar); - p->set< RCP >("QP Vector Data Layout", dl->qp_vector); - - //Output - p->set("Chemical Energy Term", "Chemical Energy Term"); - - ev = rcp(new PHAL::CahnHillChemTerm(*p)); - fm0.template registerEvaluator(ev); - } - - if(params->isParameter("Langevin Noise SD")){ - - // Form the Langevin noise term - - haveNoise = true; - - RCP p = rcp(new ParameterList("Langevin Noise Term")); - - p->set >("Parameter Library", paramLib); - - // Standard deviation of the noise - p->set("SD Value", params->get("Langevin Noise SD")); - // Time period over which to apply the noise (-1 means over the whole time) - p->set >("Langevin Noise Time Period", - params->get >("Langevin Noise Time Period", Teuchos::tuple(-1, -1))); - - //Input - p->set("Rho QP Variable Name", "Rho"); - - p->set< RCP >("QP Scalar Data Layout", dl->qp_scalar); - p->set< RCP >("QP Vector Data Layout", dl->qp_vector); - - //Output - p->set("Langevin Noise Term", "Langevin Noise Term"); - - ev = rcp(new PHAL::LangevinNoiseTerm(*p)); - fm0.template registerEvaluator(ev); - } - - { // Rho Resid - RCP p = rcp(new ParameterList("Rho Resid")); - - //Input - p->set("Weighted BF Name", "wBF"); - p->set("Weighted Gradient BF Name", "wGrad BF"); - if(haveNoise) - p->set("Langevin Noise Term", "Langevin Noise Term"); - // Accumulate in the Langevin noise term? - p->set("Have Noise", haveNoise); - - p->set("Chemical Energy Term", "Chemical Energy Term"); - p->set("Gradient QP Variable Name", "Rho Gradient"); - - p->set< RCP >("QP Scalar Data Layout", dl->qp_scalar); - p->set< RCP >("QP Vector Data Layout", dl->qp_vector); - p->set< RCP >("Node QP Scalar Data Layout", dl->node_qp_scalar); - p->set< RCP >("Node QP Vector Data Layout", dl->node_qp_vector); - - // gamma value in Equation 2.2 - p->set("gamma Value", params->get("gamma")); - - //Output - p->set("Residual Name", "Rho Residual"); - p->set< RCP >("Node Scalar Data Layout", dl->node_scalar); - - ev = rcp(new PHAL::CahnHillRhoResid(*p)); - fm0.template registerEvaluator(ev); - } - - { // W Resid - RCP p = rcp(new ParameterList("W Resid")); - - //Input - p->set("Weighted BF Name", "wBF"); - p->set("BF Name", "BF"); - p->set("Rho QP Time Derivative Variable Name", "rhoDot"); - p->set("Weighted Gradient BF Name", "wGrad BF"); - p->set("Gradient QP Variable Name", "W Gradient"); - - // Mass lump time term? - p->set("Lump Mass", params->get("Lump Mass")); - - p->set< RCP >("QP Scalar Data Layout", dl->qp_scalar); - p->set< RCP >("QP Vector Data Layout", dl->qp_vector); - p->set< RCP >("Node QP Scalar Data Layout", dl->node_qp_scalar); - p->set< RCP >("Node QP Vector Data Layout", dl->node_qp_vector); - - //Output - p->set("Residual Name", "W Residual"); - p->set< RCP >("Node Scalar Data Layout", dl->node_scalar); - - ev = rcp(new PHAL::CahnHillWResid(*p)); - fm0.template registerEvaluator(ev); - } - - if (fieldManagerChoice == Albany::BUILD_RESID_FM) { - PHX::Tag res_tag("Scatter", dl->dummy); - fm0.requireField(res_tag); - return res_tag.clone(); - } - - else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) { - Albany::ResponseUtilities respUtils(dl); - return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr); - } - - return Teuchos::null; -} - - -#endif // ALBANY_CAHNHILLPROBLEM_HPP diff --git a/src/problems/Albany_ComprNSProblem.cpp b/src/problems/Albany_ComprNSProblem.cpp deleted file mode 100644 index 9da30264e3..0000000000 --- a/src/problems/Albany_ComprNSProblem.cpp +++ /dev/null @@ -1,200 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Albany_ComprNSProblem.hpp" - -#include "Intrepid2_DefaultCubatureFactory.hpp" -#include "Shards_CellTopology.hpp" -#include "PHAL_FactoryTraits.hpp" -#include "Albany_Utils.hpp" -#include "Albany_BCUtils.hpp" -#include "Albany_ProblemUtils.hpp" -#include - - -Albany::ComprNSProblem:: -ComprNSProblem( const Teuchos::RCP& params_, - const Teuchos::RCP& paramLib_, - const int numDim_) : - Albany::AbstractProblem(params_, paramLib_), - numDim(numDim_), - use_sdbcs_(false), - params(params_) -{ - // Get number of species equations from Problem specifications - neq = params_->get("Number of PDE Equations", numDim); - -} - -Albany::ComprNSProblem:: -~ComprNSProblem() -{ -} - -void -Albany::ComprNSProblem:: -buildProblem( - Teuchos::ArrayRCP > meshSpecs, - Albany::StateManager& stateMgr) -{ - using Teuchos::rcp; - - /* Construct All Phalanx Evaluators */ - TEUCHOS_TEST_FOR_EXCEPTION(meshSpecs.size()!=1,std::logic_error,"Problem supports one Material Block"); - fm.resize(1); - fm[0] = rcp(new PHX::FieldManager); - buildEvaluators(*fm[0], *meshSpecs[0], stateMgr, BUILD_RESID_FM, - Teuchos::null); - constructDirichletEvaluators(*meshSpecs[0]); - - // Check if have Neumann sublist; throw error if attempting to specify - // Neumann BCs, but there are no sidesets in the input mesh - bool isNeumannPL = params->isSublist("Neumann BCs"); - if (isNeumannPL && !(meshSpecs[0]->ssNames.size() > 0)) { - ALBANY_ASSERT(false, "You are attempting to set Neumann BCs on a mesh with no sidesets!"); - } - - if(meshSpecs[0]->ssNames.size() > 0) // Build a sideset evaluator if sidesets are present - constructNeumannEvaluators(meshSpecs[0]); -} - -Teuchos::Array< Teuchos::RCP > -Albany::ComprNSProblem:: -buildEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList) -{ - // Call constructeEvaluators(*rfm[0], *meshSpecs[0], stateMgr); - // for each EvalT in PHAL::AlbanyTraits::BEvalTypes - ConstructEvaluatorsOp op( - *this, fm0, meshSpecs, stateMgr, fmchoice, responseList); - Sacado::mpl::for_each fe(op); - return *op.tags; -} - -void -Albany::ComprNSProblem::constructDirichletEvaluators( - const Albany::MeshSpecsStruct& meshSpecs) -{ - // Construct Dirichlet evaluators for all nodesets and names - std::vector dirichletNames(neq); - for (unsigned int i=0; i dirUtils; - dfm = dirUtils.constructBCEvaluators(meshSpecs.nsNames, dirichletNames, - this->params, this->paramLib); - use_sdbcs_ = dirUtils.useSDBCs(); - offsets_ = dirUtils.getOffsets(); - nodeSetIDs_ = dirUtils.getNodeSetIDs(); -} - -// Neumann BCs -void -Albany::ComprNSProblem::constructNeumannEvaluators(const Teuchos::RCP& meshSpecs) -{ - - std::cout << "setting Neumann evaluators" << std::endl; - // Note: we only enter this function if sidesets are defined in the mesh file - // i.e. meshSpecs.ssNames.size() > 0 - - Albany::BCUtils nbcUtils; - - // Check to make sure that Neumann BCs are given in the input file - - if(!nbcUtils.haveBCSpecified(this->params)) { - return; - } - - - // Construct BC evaluators for all side sets and names - // Note that the string index sets up the equation offset, so ordering is important - - std::vector neumannNames(neq + 1); - Teuchos::Array > offsets; - offsets.resize(neq + 1); - - neumannNames[0] = "qFluct0"; - offsets[0].resize(1); - offsets[0][0] = 0; - offsets[neq].resize(neq); - offsets[neq][0] = 0; - - if (neq>1){ - neumannNames[1] = "qFluct1"; - offsets[1].resize(1); - offsets[1][0] = 1; - offsets[neq][1] = 1; - } - - if (neq>2){ - neumannNames[2] = "qFluct2"; - offsets[2].resize(1); - offsets[2][0] = 2; - offsets[neq][2] = 2; - } - - if (neq>3){ - neumannNames[3] = "qFluct3"; - offsets[3].resize(1); - offsets[3][0] = 3; - offsets[neq][3] = 3; - } - - if (neq>4){ - neumannNames[4] = "qFluct4"; - offsets[4].resize(1); - offsets[4][0] = 4; - offsets[neq][4] = 4; - } - - neumannNames[neq] = "all"; - - // Construct BC evaluators for all possible names of conditions - // Should only specify flux vector components (dCdx, dCdy, dCdz), or dCdn, not both - std::vector condNames(2); //dCdx, dCdy, dCdz, dCdn, basal - Teuchos::ArrayRCP dof_names(1); - dof_names[0] = "qFluct"; - - // Note that sidesets are only supported for two and 3D currently - if(numDim == 2) - condNames[0] = "(dFluxdx, dFluxdy)"; - else if(numDim == 3) - condNames[0] = "(dFluxdx, dFluxdy, dFluxdz)"; - else - TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, - std::endl << "Error: Sidesets only supported in 2 and 3D." << std::endl); - - condNames[1] = "dFluxdn"; - - nfm.resize(1); // ComprNS problem only has one element block - - nfm[0] = nbcUtils.constructBCEvaluators(meshSpecs, neumannNames, dof_names, true, 0, - condNames, offsets, dl, - this->params, this->paramLib); - - -} - - -Teuchos::RCP -Albany::ComprNSProblem::getValidProblemParameters() const -{ - Teuchos::RCP validPL = - this->getGenericProblemParams("ValidComprNSProblemParams"); - - validPL->set("Number of PDE Equations", 1, "Number of PDE Equations in ComprNS equation set"); - validPL->sublist("Viscosity", false, ""); - validPL->sublist("Body Force", false, ""); - validPL->sublist("Equation Set", false, ""); - - return validPL; -} - diff --git a/src/problems/Albany_ComprNSProblem.hpp b/src/problems/Albany_ComprNSProblem.hpp deleted file mode 100644 index 10694d9516..0000000000 --- a/src/problems/Albany_ComprNSProblem.hpp +++ /dev/null @@ -1,314 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef ALBANY_COMPRNSPROBLEM_HPP -#define ALBANY_COMPRNSPROBLEM_HPP - -#include "Teuchos_RCP.hpp" -#include "Teuchos_ParameterList.hpp" - -#include "Albany_AbstractProblem.hpp" - -#include "PHAL_Workset.hpp" -#include "PHAL_Dimension.hpp" - -namespace Albany { - - /*! - * \brief Abstract interface for representing a 1-D finite element - * problem. - */ - class ComprNSProblem : public AbstractProblem { - public: - - //! Default constructor - ComprNSProblem(const Teuchos::RCP& params, - const Teuchos::RCP& paramLib, - const int numDim_); - - //! Destructor - ~ComprNSProblem(); - - //! Return number of spatial dimensions - virtual int spatialDimension() const { return numDim; } - - //! Get boolean telling code if SDBCs are utilized - virtual bool useSDBCs() const {return use_sdbcs_; } - - //! Build the PDE instantiations, boundary conditions, and initial solution - virtual void buildProblem( - Teuchos::ArrayRCP > meshSpecs, - StateManager& stateMgr); - - // Build evaluators - virtual Teuchos::Array< Teuchos::RCP > - buildEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList); - - //! Each problem must generate it's list of valid parameters - Teuchos::RCP getValidProblemParameters() const; - - private: - - //! Private to prohibit copying - ComprNSProblem(const ComprNSProblem&); - - //! Private to prohibit copying - ComprNSProblem& operator=(const ComprNSProblem&); - - public: - - //! Main problem setup routine. Not directly called, but indirectly by following functions - template Teuchos::RCP - constructEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList); - - void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs); - void constructNeumannEvaluators(const Teuchos::RCP& meshSpecs); - - protected: - int numDim; - Teuchos::RCP dl; - - /// Boolean marking whether SDBCs are used - bool use_sdbcs_; - - /// Problem params - const Teuchos::RCP params; - - }; - -} - -#include "Intrepid2_DefaultCubatureFactory.hpp" -#include "Shards_CellTopology.hpp" - -#include "Albany_Utils.hpp" -#include "Albany_ProblemUtils.hpp" -#include "Albany_EvaluatorUtils.hpp" -#include "Albany_ResponseUtilities.hpp" - -#include "PHAL_DOFVecGradInterpolation.hpp" - -#include "PHAL_ComprNSResid.hpp" - -#include "PHAL_Neumann.hpp" -#include "PHAL_Source.hpp" -#include "PHAL_ComprNSBodyForce.hpp" -#include "PHAL_ComprNSViscosity.hpp" - -template -Teuchos::RCP -Albany::ComprNSProblem::constructEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fieldManagerChoice, - const Teuchos::RCP& responseList) -{ - using Teuchos::RCP; - using Teuchos::rcp; - using Teuchos::ParameterList; - using PHX::DataLayout; - using PHX::MDALayout; - using std::vector; - using std::string; - using std::map; - using PHAL::AlbanyTraits; - - RCP > - intrepidBasis = Albany::getIntrepid2Basis(meshSpecs.ctd); - RCP cellType = rcp(new shards::CellTopology (&meshSpecs.ctd)); - - // Problem is transient - TEUCHOS_TEST_FOR_EXCEPTION( - number_of_time_deriv != 1, - std::logic_error, - "Albany_ComprNSProblem must be defined as a transient calculation."); - - const int numNodes = intrepidBasis->getCardinality(); - const int worksetSize = meshSpecs.worksetSize; - - Intrepid2::DefaultCubatureFactory cubFactory; - RCP > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree); - - const int numQPts = cubature->getNumPoints(); - const int numVertices = cellType->getNodeCount(); - - *out << "Field Dimensions: Workset=" << worksetSize - << ", Vertices= " << numVertices - << ", Nodes= " << numNodes - << ", QuadPts= " << numQPts - << ", Dim= " << numDim << std::endl; - - int vecDim = neq; - - dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim, vecDim)); - Albany::EvaluatorUtils evalUtils(dl); - int offset=0; - - // Temporary variable used numerous times below - Teuchos::RCP > ev; - - // Define Field Names - - Teuchos::ArrayRCP dof_names(1); - Teuchos::ArrayRCP dof_names_dot(1); - Teuchos::ArrayRCP resid_names(1); - dof_names[0] = "qFluct"; - dof_names_dot[0] = dof_names[0]+"_dot"; - resid_names[0] = "ComprNS Residual"; - fm0.template registerEvaluator - (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dot, offset)); - - fm0.template registerEvaluator - (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], offset)); - - fm0.template registerEvaluator - (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dot[0], offset)); - - // fm0.template registerEvaluator - // (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], offset)); - - fm0.template registerEvaluator - (evalUtils.constructScatterResidualEvaluator(true, resid_names,offset, "Scatter ComprNS")); - - fm0.template registerEvaluator - (evalUtils.constructGatherCoordinateVectorEvaluator()); - - fm0.template registerEvaluator - (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature)); - - fm0.template registerEvaluator - (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature)); - - { // Specialized DofVecGrad Interpolation for this problem - - RCP p = rcp(new ParameterList("DOFVecGrad Interpolation "+dof_names[0])); - // Input - p->set("Variable Name", dof_names[0]); - - p->set("Gradient BF Name", "Grad BF"); - p->set("Offset of First DOF", offset); - - // Output (assumes same Name as input) - p->set("Gradient Variable Name", dof_names[0]+" Gradient"); - - ev = rcp(new PHAL::DOFVecGradInterpolation(*p,dl)); - fm0.template registerEvaluator(ev); - } - - { // ComprNS Resid - RCP p = rcp(new ParameterList("ComprNS Resid")); - - //Input - p->set("Weighted BF Name", "wBF"); - p->set("Weighted Gradient BF Name", "wGrad BF"); - p->set("QP Variable Name", "qFluct"); - p->set("QP Time Derivative Variable Name", "qFluct_dot"); - p->set("Gradient QP Variable Name", "qFluct Gradient"); - p->set("Body Force Name", "Body Force"); - - p->set< RCP >("QP Vector Data Layout", dl->qp_vector); - p->set< RCP >("QP Tensor Data Layout", dl->qp_vecgradient); - p->set< RCP >("Node QP Scalar Data Layout", dl->node_qp_scalar); - p->set< RCP >("Node QP Gradient Data Layout", dl->node_qp_gradient); - - p->set("Viscosity Mu QP Variable Name", "Viscosity Mu"); - p->set("Viscosity Lambda QP Variable Name", "Viscosity Lambda"); - p->set("Viscosity Kappa QP Variable Name", "Viscosity Kappa"); - p->set< RCP >("QP Scalar Data Layout", dl->qp_scalar); - - p->set("Stress Tau11 QP Variable Name", "Stress Tau11"); - p->set("Stress Tau12 QP Variable Name", "Stress Tau12"); - p->set("Stress Tau13 QP Variable Name", "Stress Tau13"); - p->set("Stress Tau22 QP Variable Name", "Stress Tau22"); - p->set("Stress Tau23 QP Variable Name", "Stress Tau23"); - p->set("Stress Tau33 QP Variable Name", "Stress Tau33"); - - p->set >("Parameter Library", paramLib); - Teuchos::ParameterList& paramList = params->sublist("Equation Set"); - p->set("Parameter List", ¶mList); - - //Output - p->set("Residual Name", "ComprNS Residual"); - p->set< RCP >("Node Vector Data Layout", dl->node_vector); - - ev = rcp(new PHAL::ComprNSResid(*p)); - fm0.template registerEvaluator(ev); - } - { // Viscosity - RCP p = rcp(new ParameterList("Viscosity")); - - //Input - p->set("Coordinate Vector Name", "Coord Vec"); - p->set< RCP >("QP Vector Data Layout", dl->qp_vector); - p->set< RCP >("QP Gradient Data Layout", dl->qp_gradient); - p->set< RCP >("QP Tensor Data Layout", dl->qp_vecgradient); - p->set("QP Variable Name", "qFluct"); - p->set("Gradient QP Variable Name", "qFluct Gradient"); - - p->set >("Parameter Library", paramLib); - Teuchos::ParameterList& paramList = params->sublist("Viscosity"); - p->set("Parameter List", ¶mList); - - //Output - p->set("Viscosity Mu QP Variable Name", "Viscosity Mu"); - p->set("Viscosity Lambda QP Variable Name", "Viscosity Lambda"); - p->set("Viscosity Kappa QP Variable Name", "Viscosity Kappa"); - p->set< RCP >("QP Scalar Data Layout", dl->qp_scalar); - p->set("Stress Tau11 QP Variable Name", "Stress Tau11"); - p->set("Stress Tau12 QP Variable Name", "Stress Tau12"); - p->set("Stress Tau13 QP Variable Name", "Stress Tau13"); - p->set("Stress Tau22 QP Variable Name", "Stress Tau22"); - p->set("Stress Tau23 QP Variable Name", "Stress Tau23"); - p->set("Stress Tau33 QP Variable Name", "Stress Tau33"); - - ev = rcp(new PHAL::ComprNSViscosity(*p)); - fm0.template registerEvaluator(ev); - - } - - { // ComprNS Body Force - RCP p = rcp(new ParameterList("Body Force")); - - //Input - p->set< RCP >("QP Scalar Data Layout", dl->qp_scalar); - p->set< RCP >("QP Vector Data Layout", dl->qp_vector); - p->set< RCP >("QP Gradient Data Layout", dl->qp_gradient); - p->set("Coordinate Vector Name", "Coord Vec"); - - Teuchos::ParameterList& paramList = params->sublist("Body Force"); - p->set("Parameter List", ¶mList); - - //Output - p->set("Body Force Name", "Body Force"); - - ev = rcp(new PHAL::ComprNSBodyForce(*p)); - fm0.template registerEvaluator(ev); - } - - - if (fieldManagerChoice == Albany::BUILD_RESID_FM) { - PHX::Tag res_tag("Scatter ComprNS", dl->dummy); - fm0.requireField(res_tag); - } - else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) { - Albany::ResponseUtilities respUtils(dl); - return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr); - } - - return Teuchos::null; -} -#endif // ALBANY_COMPRNSPROBLEM_HPP diff --git a/src/problems/Albany_DemoProblemFactory.cpp b/src/problems/Albany_DemoProblemFactory.cpp index 1ac83f852a..a92c6c027e 100644 --- a/src/problems/Albany_DemoProblemFactory.cpp +++ b/src/problems/Albany_DemoProblemFactory.cpp @@ -6,15 +6,11 @@ #include "Albany_DemoProblemFactory.hpp" -#include "Albany_CahnHillProblem.hpp" #include "Albany_Helmholtz2DProblem.hpp" #include "Albany_NavierStokes.hpp" -#include "Albany_LinComprNSProblem.hpp" #include "Albany_AdvDiffProblem.hpp" #include "Albany_ReactDiffSystem.hpp" -#include "Albany_ComprNSProblem.hpp" #include "Albany_ODEProblem.hpp" -#include "Albany_PNPProblem.hpp" #include "Albany_ThermoElectrostaticsProblem.hpp" #include "Albany_ThermalProblem.hpp" #include "Albany_AdvectionProblem.hpp" @@ -37,24 +33,15 @@ int getNumDim(std::string const& key) bool DemoProblemFactory::provides (const std::string& key) const { - return key == "CahnHill 2D" || - key == "ODE" || + return key == "ODE" || key == "Helmholtz 2D" || key == "NavierStokes 1D" || key == "NavierStokes 2D" || key == "NavierStokes 3D" || - key == "LinComprNS 1D" || key == "AdvDiff 1D" || key == "AdvDiff 2D" || key=="Reaction-Diffusion System 3D" || key == "Reaction-Diffusion System" || - key == "LinComprNS 2D" || - key == "LinComprNS 3D" || - key == "ComprNS 2D" || - key == "ComprNS 3D" || - key == "PNP 1D" || - key == "PNP 2D" || - key == "PNP 3D" || key == "Thermal 1D" || key == "Thermal 2D" || key == "Thermal 3D" || @@ -78,9 +65,7 @@ create (const std::string& key, auto problemParams = Teuchos::sublist(topLevelParams, "Problem", true); auto discretizationParams = Teuchos::sublist(topLevelParams, "Discretization"); - if (key == "CahnHill 2D") { - problem = Teuchos::rcp(new CahnHillProblem(problemParams, paramLib, 2, comm)); - } else if (key == "ODE") { + if (key == "ODE") { problem = Teuchos::rcp(new ODEProblem(problemParams, paramLib, 0)); } else if (key == "Helmholtz 2D") { problem = Teuchos::rcp(new Helmholtz2DProblem(problemParams, paramLib)); @@ -90,8 +75,6 @@ create (const std::string& key, problem = Teuchos::rcp(new NavierStokes(problemParams, paramLib, 2)); } else if (key == "NavierStokes 3D") { problem = Teuchos::rcp(new NavierStokes(problemParams, paramLib, 3)); - } else if (key == "LinComprNS 1D") { - problem = Teuchos::rcp(new LinComprNSProblem(problemParams, paramLib, 1)); } else if (key == "AdvDiff 1D") { problem = Teuchos::rcp(new AdvDiffProblem(problemParams, paramLib, 1)); } else if (key == "AdvDiff 2D") { @@ -99,20 +82,6 @@ create (const std::string& key, } else if (key=="Reaction-Diffusion System 3D" || key == "Reaction-Diffusion System") { problem = Teuchos::rcp(new ReactDiffSystem(problemParams, paramLib, 3)); - } else if (key == "LinComprNS 2D") { - problem = Teuchos::rcp(new LinComprNSProblem(problemParams, paramLib, 2)); - } else if (key == "LinComprNS 3D") { - problem = Teuchos::rcp(new LinComprNSProblem(problemParams, paramLib, 3)); - } else if (key == "ComprNS 2D") { - problem = Teuchos::rcp(new ComprNSProblem(problemParams, paramLib, 2)); - } else if (key == "ComprNS 3D") { - problem = Teuchos::rcp(new ComprNSProblem(problemParams, paramLib, 3)); - } else if (key == "PNP 1D") { - problem = Teuchos::rcp(new PNPProblem(problemParams, paramLib, 1)); - } else if (key == "PNP 2D") { - problem = Teuchos::rcp(new PNPProblem(problemParams, paramLib, 2)); - } else if (key == "PNP 3D") { - problem = Teuchos::rcp(new PNPProblem(problemParams, paramLib, 3)); } else if (key == "ThermoElectrostatics 1D") { problem = Teuchos::rcp(new ThermoElectrostaticsProblem(problemParams, paramLib, 1)); } else if (key == "ThermoElectrostatics 2D") { diff --git a/src/problems/Albany_LinComprNSProblem.cpp b/src/problems/Albany_LinComprNSProblem.cpp deleted file mode 100644 index 11dc8fddd5..0000000000 --- a/src/problems/Albany_LinComprNSProblem.cpp +++ /dev/null @@ -1,99 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Albany_LinComprNSProblem.hpp" - -#include "Intrepid2_DefaultCubatureFactory.hpp" -#include "Shards_CellTopology.hpp" -#include "PHAL_FactoryTraits.hpp" -#include "Albany_Utils.hpp" -#include "Albany_BCUtils.hpp" -#include "Albany_ProblemUtils.hpp" -#include - - -Albany::LinComprNSProblem:: -LinComprNSProblem( const Teuchos::RCP& params_, - const Teuchos::RCP& paramLib_, - const int numDim_) : - Albany::AbstractProblem(params_, paramLib_), - numDim(numDim_), - use_sdbcs_(false) -{ - // Get number of species equations from Problem specifications - neq = params_->get("Number of PDE Equations", numDim); -} - -Albany::LinComprNSProblem:: -~LinComprNSProblem() -{ -} - -void -Albany::LinComprNSProblem:: -buildProblem( - Teuchos::ArrayRCP > meshSpecs, - Albany::StateManager& stateMgr) -{ - using Teuchos::rcp; - - /* Construct All Phalanx Evaluators */ - TEUCHOS_TEST_FOR_EXCEPTION(meshSpecs.size()!=1,std::logic_error,"Problem supports one Material Block"); - fm.resize(1); - fm[0] = rcp(new PHX::FieldManager); - buildEvaluators(*fm[0], *meshSpecs[0], stateMgr, BUILD_RESID_FM, - Teuchos::null); - constructDirichletEvaluators(*meshSpecs[0]); -} - -Teuchos::Array< Teuchos::RCP > -Albany::LinComprNSProblem:: -buildEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList) -{ - // Call constructeEvaluators(*rfm[0], *meshSpecs[0], stateMgr); - // for each EvalT in PHAL::AlbanyTraits::BEvalTypes - ConstructEvaluatorsOp op( - *this, fm0, meshSpecs, stateMgr, fmchoice, responseList); - Sacado::mpl::for_each fe(op); - return *op.tags; -} - -void -Albany::LinComprNSProblem::constructDirichletEvaluators( - const Albany::MeshSpecsStruct& meshSpecs) -{ - // Construct Dirichlet evaluators for all nodesets and names - std::vector dirichletNames(neq); - for (unsigned int i=0; i dirUtils; - dfm = dirUtils.constructBCEvaluators(meshSpecs.nsNames, dirichletNames, - this->params, this->paramLib); - use_sdbcs_ = dirUtils.useSDBCs(); - offsets_ = dirUtils.getOffsets(); - nodeSetIDs_ = dirUtils.getNodeSetIDs(); -} - -Teuchos::RCP -Albany::LinComprNSProblem::getValidProblemParameters() const -{ - Teuchos::RCP validPL = - this->getGenericProblemParams("ValidLinComprNSProblemParams"); - - validPL->set("Number of PDE Equations", 1, "Number of PDE Equations in LinComprNS equation set"); - validPL->sublist("Body Force", false, ""); - validPL->sublist("Equation Set", false, ""); - - return validPL; -} - diff --git a/src/problems/Albany_LinComprNSProblem.hpp b/src/problems/Albany_LinComprNSProblem.hpp deleted file mode 100644 index c0e3cddef2..0000000000 --- a/src/problems/Albany_LinComprNSProblem.hpp +++ /dev/null @@ -1,272 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef ALBANY_LINCOMPRNSPROBLEM_HPP -#define ALBANY_LINCOMPRNSPROBLEM_HPP - -#include "Teuchos_RCP.hpp" -#include "Teuchos_ParameterList.hpp" - -#include "Albany_AbstractProblem.hpp" - -#include "PHAL_Workset.hpp" -#include "PHAL_Dimension.hpp" - -namespace Albany { - - /*! - * \brief Abstract interface for representing a 1-D finite element - * problem. - */ - class LinComprNSProblem : public AbstractProblem { - public: - - //! Default constructor - LinComprNSProblem(const Teuchos::RCP& params, - const Teuchos::RCP& paramLib, - const int numDim_); - - //! Destructor - ~LinComprNSProblem(); - - //! Return number of spatial dimensions - virtual int spatialDimension() const { return numDim; } - - //! Get boolean telling code if SDBCs are utilized - virtual bool useSDBCs() const {return use_sdbcs_; } - - //! Build the PDE instantiations, boundary conditions, and initial solution - virtual void buildProblem( - Teuchos::ArrayRCP > meshSpecs, - StateManager& stateMgr); - - // Build evaluators - virtual Teuchos::Array< Teuchos::RCP > - buildEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList); - - //! Each problem must generate it's list of valid parameters - Teuchos::RCP getValidProblemParameters() const; - - private: - - //! Private to prohibit copying - LinComprNSProblem(const LinComprNSProblem&); - - //! Private to prohibit copying - LinComprNSProblem& operator=(const LinComprNSProblem&); - - public: - - //! Main problem setup routine. Not directly called, but indirectly by following functions - template Teuchos::RCP - constructEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList); - - void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs); - - protected: - int numDim; - - /// Boolean marking whether SDBCs are used - bool use_sdbcs_; - - }; - -} - -#include "Intrepid2_DefaultCubatureFactory.hpp" -#include "Shards_CellTopology.hpp" - -#include "Albany_Utils.hpp" -#include "Albany_ProblemUtils.hpp" -#include "Albany_EvaluatorUtils.hpp" -#include "Albany_ResponseUtilities.hpp" - -#include "PHAL_DOFVecGradInterpolation.hpp" - -#include "PHAL_LinComprNSResid.hpp" - -#include "PHAL_Source.hpp" -#include "PHAL_LinComprNSBodyForce.hpp" - -template -Teuchos::RCP -Albany::LinComprNSProblem::constructEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fieldManagerChoice, - const Teuchos::RCP& responseList) -{ - using Teuchos::RCP; - using Teuchos::rcp; - using Teuchos::ParameterList; - using PHX::DataLayout; - using PHX::MDALayout; - using std::vector; - using std::string; - using std::map; - using PHAL::AlbanyTraits; - - RCP > - intrepidBasis = Albany::getIntrepid2Basis(meshSpecs.ctd); - RCP cellType = rcp(new shards::CellTopology (&meshSpecs.ctd)); - - const int numNodes = intrepidBasis->getCardinality(); - const int worksetSize = meshSpecs.worksetSize; - - Intrepid2::DefaultCubatureFactory cubFactory; - RCP > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree); - - const int numQPts = cubature->getNumPoints(); - const int numVertices = cellType->getNodeCount(); - - *out << "Field Dimensions: Workset=" << worksetSize - << ", Vertices= " << numVertices - << ", Nodes= " << numNodes - << ", QuadPts= " << numQPts - << ", Dim= " << numDim << std::endl; - - int vecDim = neq; - - RCP dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim, vecDim)); - Albany::EvaluatorUtils evalUtils(dl); - int offset=0; - - // Make sure we are transient - TEUCHOS_TEST_FOR_EXCEPTION( - number_of_time_deriv < 0 || number_of_time_deriv > 1, - std::logic_error, - "Albany_LinComprNSProblem must be defined as a steady or transient calculation."); - - // Temporary variable used numerous times below - Teuchos::RCP > ev; - - // Define Field Names - - Teuchos::ArrayRCP dof_names(1); - Teuchos::ArrayRCP dof_names_dot(1); - Teuchos::ArrayRCP resid_names(1); - dof_names[0] = "qFluct"; - if(number_of_time_deriv == 1) - dof_names_dot[0] = dof_names[0]+"_dot"; - resid_names[0] = "LinComprNS Residual"; - if(number_of_time_deriv == 1) - fm0.template registerEvaluator - (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dot, offset)); - else - fm0.template registerEvaluator - (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names, offset)); - - - fm0.template registerEvaluator - (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], offset)); - - if(number_of_time_deriv == 1) - fm0.template registerEvaluator - (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dot[0], offset)); - - // fm0.template registerEvaluator - // (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], offset)); - - fm0.template registerEvaluator - (evalUtils.constructScatterResidualEvaluator(true, resid_names, offset, "Scatter LinComprNS")); - - fm0.template registerEvaluator - (evalUtils.constructGatherCoordinateVectorEvaluator()); - - fm0.template registerEvaluator - (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature)); - - fm0.template registerEvaluator - (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature)); - - { // Specialized DofVecGrad Interpolation for this problem - - RCP p = rcp(new ParameterList("DOFVecGrad Interpolation "+dof_names[0])); - // Input - p->set("Variable Name", dof_names[0]); - p->set("Gradient BF Name", "Grad BF"); - p->set("Offset of First DOF", offset); - - // Output (assumes same Name as input) - p->set("Gradient Variable Name", dof_names[0]+" Gradient"); - - ev = rcp(new PHAL::DOFVecGradInterpolation(*p,dl)); - fm0.template registerEvaluator(ev); - } - - { // LinComprNS Resid - RCP p = rcp(new ParameterList("LinComprNS Resid")); - - //Input - p->set("Weighted BF Name", "wBF"); - p->set("Weighted Gradient BF Name", "wGrad BF"); - p->set("QP Variable Name", "qFluct"); - p->set("QP Time Derivative Variable Name", "qFluct_dot"); - p->set("Gradient QP Variable Name", "qFluct Gradient"); - p->set("Body Force Name", "Body Force"); - if(number_of_time_deriv == 0) - p->set("Disable Transient", true); - - p->set< RCP >("QP Vector Data Layout", dl->qp_vector); - p->set< RCP >("QP Tensor Data Layout", dl->qp_vecgradient); - p->set< RCP >("Node QP Scalar Data Layout", dl->node_qp_scalar); - p->set< RCP >("Node QP Gradient Data Layout", dl->node_qp_gradient); - - p->set >("Parameter Library", paramLib); - Teuchos::ParameterList& paramList = params->sublist("Equation Set"); - p->set("Parameter List", ¶mList); - - //Output - p->set("Residual Name", "LinComprNS Residual"); - p->set< RCP >("Node Vector Data Layout", dl->node_vector); - - ev = rcp(new PHAL::LinComprNSResid(*p)); - fm0.template registerEvaluator(ev); - } - - { // LinComprNS Body Force - RCP p = rcp(new ParameterList("Body Force")); - - //Input - p->set< RCP >("QP Scalar Data Layout", dl->qp_scalar); - p->set< RCP >("QP Vector Data Layout", dl->qp_vector); - p->set< RCP >("QP Gradient Data Layout", dl->qp_gradient); - p->set("Coordinate Vector Name", "Coord Vec"); - - Teuchos::ParameterList& paramList = params->sublist("Body Force"); - p->set("Parameter List", ¶mList); - - //Output - p->set("Body Force Name", "Body Force"); - - ev = rcp(new PHAL::LinComprNSBodyForce(*p)); - fm0.template registerEvaluator(ev); - } - - - if (fieldManagerChoice == Albany::BUILD_RESID_FM) { - PHX::Tag res_tag("Scatter LinComprNS", dl->dummy); - fm0.requireField(res_tag); - } - else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) { - Albany::ResponseUtilities respUtils(dl); - return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr); - } - - return Teuchos::null; -} -#endif // ALBANY_LinComprNSPROBLEM_HPP diff --git a/src/problems/Albany_PNPProblem.cpp b/src/problems/Albany_PNPProblem.cpp deleted file mode 100644 index c0ec4d2642..0000000000 --- a/src/problems/Albany_PNPProblem.cpp +++ /dev/null @@ -1,178 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Albany_PNPProblem.hpp" - -#include "Intrepid2_DefaultCubatureFactory.hpp" -#include "Shards_CellTopology.hpp" -#include "PHAL_FactoryTraits.hpp" -#include "Albany_Utils.hpp" -#include "Albany_BCUtils.hpp" -#include "Albany_ProblemUtils.hpp" - -Albany::PNPProblem:: -PNPProblem( const Teuchos::RCP& params_, - const Teuchos::RCP& paramLib_, - const int numDim_) : - Albany::AbstractProblem(params_, paramLib_), - numDim(numDim_), - use_sdbcs_(false), - params(params_) -{ - - // Compute number of equations - numSpecies = params->get("Number of Species", 1); - int num_eq = numSpecies + 1; - this->setNumEquations(num_eq); - - // Print out a summary of the problem - *out << "PNP problem: with numSpecies = " << numSpecies << std::endl; - -} - -Albany::PNPProblem:: -~PNPProblem() -{ -} - -void -Albany::PNPProblem:: -buildProblem( - Teuchos::ArrayRCP > meshSpecs, - Albany::StateManager& stateMgr) -{ - using Teuchos::rcp; - /* Construct All Phalanx Evaluators */ - int physSets = meshSpecs.size(); - std::cout << "PNP Problem Num MeshSpecs: " << physSets << std::endl; - fm.resize(physSets); - - for (int ps=0; ps); - buildEvaluators(*fm[ps], *meshSpecs[ps], stateMgr, BUILD_RESID_FM, - Teuchos::null); - } - - if(meshSpecs[0]->nsNames.size() > 0) // Build a nodeset evaluator if nodesets are present - constructDirichletEvaluators(*meshSpecs[0]); - - - // Check if have Neumann sublist; throw error if attempting to specify - // Neumann BCs, but there are no sidesets in the input mesh - bool isNeumannPL = params->isSublist("Neumann BCs"); - if (isNeumannPL && !(meshSpecs[0]->ssNames.size() > 0)) { - ALBANY_ASSERT(false, "You are attempting to set Neumann BCs on a mesh with no sidesets!"); - } - - - if(meshSpecs[0]->ssNames.size() > 0) // Build a sideset evaluator if sidesets are present - constructNeumannEvaluators(meshSpecs[0]); -} - -Teuchos::Array< Teuchos::RCP > -Albany::PNPProblem:: -buildEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList) -{ - // Call constructeEvaluators(*rfm[0], *meshSpecs[0], stateMgr); - // for each EvalT in PHAL::AlbanyTraits::BEvalTypes - Albany::ConstructEvaluatorsOp op( - *this, fm0, meshSpecs, stateMgr, fmchoice, responseList); - Sacado::mpl::for_each fe(op); - return *op.tags; -} - -void -Albany::PNPProblem::constructDirichletEvaluators( - const Albany::MeshSpecsStruct& meshSpecs) -{ - // Construct Dirichlet evaluators for all nodesets and names - std::vector dirichletNames(neq); - int idx = 0; - for(idx = 0; idx dirUtils; - dfm = dirUtils.constructBCEvaluators(meshSpecs.nsNames, dirichletNames, - this->params, this->paramLib); - use_sdbcs_ = dirUtils.useSDBCs(); - offsets_ = dirUtils.getOffsets(); - nodeSetIDs_ = dirUtils.getNodeSetIDs(); -} - -//Neumann BCs -void -Albany::PNPProblem::constructNeumannEvaluators(const Teuchos::RCP& meshSpecs) -{ - - // Note: we only enter this function if sidesets are defined in the mesh file - // i.e. meshSpecs.ssNames.size() > 0 - - Albany::BCUtils nbcUtils; - - // Check to make sure that Neumann BCs are given in the input file - - if(!nbcUtils.haveBCSpecified(this->params)) { - return; - } - - // Construct BC evaluators for all side sets and names - // Note that the string index sets up the equation offset, so ordering is important - // - // Currently we aren't exactly doing this right. I think to do this - // correctly we need different neumann evaluators for each DOF (velocity, - // pressure, temperature, flux) since velocity is a vector and the - // others are scalars. The dof_names stuff is only used - // for robin conditions, so at this point, as long as we don't enable - // robin conditions, this should work. - - std::vector nbcNames; - Teuchos::RCP< Teuchos::Array > dof_names = - Teuchos::rcp(new Teuchos::Array); - Teuchos::Array > offsets; - int idx = 0; - for(idx = 0; idx(1,idx)); - } - nbcNames.push_back("Phi"); - offsets.push_back(Teuchos::Array(1,idx++)); - dof_names->push_back("Concentration"); - dof_names->push_back("Potential"); - - // Construct BC evaluators for all possible names of conditions - // Should only specify flux vector components (dudx, dudy, dudz), or dudn, not both - std::vector condNames; //dudx, dudy, dudz, dudn, basal - - - nfm[0] = nbcUtils.constructBCEvaluators(meshSpecs, nbcNames, - Teuchos::arcp(dof_names), - true, 0, condNames, offsets, dl, - this->params, this->paramLib); -} - - - -Teuchos::RCP -Albany::PNPProblem::getValidProblemParameters() const -{ - Teuchos::RCP validPL = - this->getGenericProblemParams("ValidPNPParams"); - - validPL->sublist("Body Force", false, ""); - validPL->sublist("Permittivity", false, ""); - validPL->set("Number of Species", 1, "Number of diffusing species"); - - return validPL; -} - diff --git a/src/problems/Albany_PNPProblem.hpp b/src/problems/Albany_PNPProblem.hpp deleted file mode 100644 index 07efbf8a82..0000000000 --- a/src/problems/Albany_PNPProblem.hpp +++ /dev/null @@ -1,304 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef Albany_PNPPROBLEM_HPP -#define Albany_PNPPROBLEM_HPP - -#include "Teuchos_RCP.hpp" -#include "Teuchos_ParameterList.hpp" - -#include "Albany_AbstractProblem.hpp" - -#include "PHAL_Workset.hpp" -#include "PHAL_Dimension.hpp" - -namespace Albany { - - /*! - * \brief Abstract interface for representing a 1-D finite element - * problem. - */ - class PNPProblem : public Albany::AbstractProblem { - public: - - //! Default constructor - PNPProblem(const Teuchos::RCP& params, - const Teuchos::RCP& paramLib, - const int numDim_); - - //! Destructor - ~PNPProblem(); - - //! Build the PDE instantiations, boundary conditions, and initial solution - virtual void buildProblem( - Teuchos::ArrayRCP > meshSpecs, - Albany::StateManager& stateMgr); - - //! Return number of spatial dimensions - virtual int spatialDimension() const { return numDim; } - - //! Get boolean telling code if SDBCs are utilized - virtual bool useSDBCs() const {return use_sdbcs_; } - - // Build evaluators - virtual Teuchos::Array< Teuchos::RCP > - buildEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList); - - //! Each problem must generate it's list of valid parameters - Teuchos::RCP getValidProblemParameters() const; - - private: - - //! Private to prohibit copying - PNPProblem(const PNPProblem&); - - //! Private to prohibit copying - PNPProblem& operator=(const PNPProblem&); - - public: - - //! Main problem setup routine. Not directly called, but indirectly by following functions - template Teuchos::RCP - constructEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fmchoice, - const Teuchos::RCP& responseList); - - void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs); - void constructNeumannEvaluators(const Teuchos::RCP& meshSpecs); - - - protected: - - int numDim; //! number of spatial dimensions - int numSpecies; //! number of species - - Teuchos::RCP dl; - - /// Boolean marking whether SDBCs are used - bool use_sdbcs_; - - //! Problem PL - const Teuchos::RCP params; - - - }; - -} - -#include "Intrepid2_DefaultCubatureFactory.hpp" -#include "Shards_CellTopology.hpp" - -#include "Albany_Utils.hpp" -#include "Albany_ProblemUtils.hpp" -#include "Albany_EvaluatorUtils.hpp" -#include "Albany_ResponseUtilities.hpp" - -#include "PNP_PotentialResid.hpp" -#include "PNP_ConcentrationResid.hpp" -#include "PHAL_Permittivity.hpp" -#include "PHAL_Neumann.hpp" - -template -Teuchos::RCP -Albany::PNPProblem::constructEvaluators( - PHX::FieldManager& fm0, - const Albany::MeshSpecsStruct& meshSpecs, - Albany::StateManager& stateMgr, - Albany::FieldManagerChoice fieldManagerChoice, - const Teuchos::RCP& responseList) -{ - using Teuchos::RCP; - using Teuchos::rcp; - using Teuchos::ParameterList; - using PHX::DataLayout; - using PHX::MDALayout; - using std::vector; - using std::string; - using std::map; - using PHAL::AlbanyTraits; - - RCP > - intrepidBasis = Albany::getIntrepid2Basis(meshSpecs.ctd); - RCP cellType = rcp(new shards::CellTopology (&meshSpecs.ctd)); - - const int numNodes = intrepidBasis->getCardinality(); - const int worksetSize = meshSpecs.worksetSize; - - Intrepid2::DefaultCubatureFactory cubFactory; - RCP > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree); - - const int numQPts = cubature->getNumPoints(); - const int numVertices = cellType->getNodeCount(); - - *out << "Field Dimensions: Workset=" << worksetSize - << ", Vertices= " << numVertices - << ", Nodes= " << numNodes - << ", QuadPts= " << numQPts - << ", Dim= " << numDim << std::endl; - - - dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim, numSpecies)); - Albany::EvaluatorUtils evalUtils(dl); - int offset=0; - - // Problem is transient - TEUCHOS_TEST_FOR_EXCEPTION( - number_of_time_deriv != 0, -// number_of_time_deriv < 0 || number_of_time_deriv > 1, - std::logic_error, - "Albany_PNPProblem must be defined as a steady calculation."); -// "Albany_PNPProblem must be defined as a steady or transient calculation."); - - // Temporary variable used numerous times below - Teuchos::RCP > ev; - - // Define Field Names - - - { - Teuchos::ArrayRCP dof_names(1); - Teuchos::ArrayRCP dof_names_dot(1); - Teuchos::ArrayRCP resid_names(1); - dof_names[0] = "Potential"; - resid_names[0] = "Potential Residual"; - if(number_of_time_deriv > 0) - fm0.template registerEvaluator - (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot, offset)); - else - fm0.template registerEvaluator - (evalUtils.constructGatherSolutionEvaluator_noTransient(false, dof_names, offset)); - - fm0.template registerEvaluator - (evalUtils.constructDOFInterpolationEvaluator(dof_names[0], offset)); - - fm0.template registerEvaluator - (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[0], offset)); - - fm0.template registerEvaluator - (evalUtils.constructScatterResidualEvaluator(false, resid_names,offset, "Scatter Potential")); - offset ++; - } - - { - Teuchos::ArrayRCP dof_names(1); - Teuchos::ArrayRCP dof_names_dot(1); - Teuchos::ArrayRCP resid_names(1); - dof_names[0] = "Concentration"; - if(number_of_time_deriv > 0) - dof_names_dot[0] = dof_names[0]+"_dot"; - resid_names[0] = "Concentration Residual"; - if(number_of_time_deriv > 0) - fm0.template registerEvaluator - (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dot, offset)); - else - fm0.template registerEvaluator - (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names, offset)); - - fm0.template registerEvaluator - (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], offset)); - - if(number_of_time_deriv > 0) - fm0.template registerEvaluator - (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dot[0], offset)); - - fm0.template registerEvaluator - (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], offset)); - - fm0.template registerEvaluator - (evalUtils.constructScatterResidualEvaluator(true, resid_names,offset, "Scatter Concentration")); - offset += numSpecies; - } - - fm0.template registerEvaluator - (evalUtils.constructGatherCoordinateVectorEvaluator()); - - fm0.template registerEvaluator - (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature)); - - fm0.template registerEvaluator - (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature)); - - { // Permittivity - RCP p = rcp(new ParameterList); - - p->set("QP Variable Name", "Permittivity"); - p->set("QP Coordinate Vector Name", "Coord Vec"); - p->set< RCP >("Node Data Layout", dl->node_scalar); - p->set< RCP >("QP Scalar Data Layout", dl->qp_scalar); - p->set< RCP >("QP Vector Data Layout", dl->qp_vector); - - p->set >("Parameter Library", paramLib); - Teuchos::ParameterList& paramList = params->sublist("Permittivity"); - p->set("Parameter List", ¶mList); - - ev = rcp(new PHAL::Permittivity(*p)); - fm0.template registerEvaluator(ev); - } - - { // Concentration Resid - RCP p = rcp(new ParameterList("Concentration Resid")); - - //Input - p->set("Weighted BF Name", "wBF"); - p->set("Weighted Gradient BF Name", "wGrad BF"); - if(number_of_time_deriv == 0) - p->set("Disable Transient", true); - - // Variable names hardwired to Concentration and Potential in evaluators - - p->set >("Parameter Library", paramLib); - - //Output - p->set("Residual Name", "Concentration Residual"); - - ev = rcp(new PNP::ConcentrationResid(*p,dl)); - fm0.template registerEvaluator(ev); - } - - - { // Potential Resid - RCP p = rcp(new ParameterList("Potential Resid")); - - //Input - p->set("Weighted BF Name", "wBF"); - p->set("Permittivity Name", "Permittivity"); - p->set< RCP >("QP Scalar Data Layout", dl->qp_scalar); - p->set("Weighted Gradient BF Name", "wGrad BF"); - // Variable names hardwired to Concentration and Potential in evaluators - - //Output - p->set("Residual Name", "Potential Residual"); - - ev = rcp(new PNP::PotentialResid(*p,dl)); - fm0.template registerEvaluator(ev); - } - - - if (fieldManagerChoice == Albany::BUILD_RESID_FM) { - Teuchos::RCP ret_tag; - PHX::Tag con_tag("Scatter Potential", dl->dummy); - fm0.requireField(con_tag); - PHX::Tag mom_tag("Scatter Concentration", dl->dummy); - fm0.requireField(mom_tag); - ret_tag = mom_tag.clone(); - return ret_tag; - } - else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) { - Albany::ResponseUtilities respUtils(dl); - return respUtils.constructResponses(fm0, *responseList, stateMgr); - } - - return Teuchos::null; -} -#endif // Albany_PNP_HPP diff --git a/src/problems/Albany_ProblemUtils.cpp b/src/problems/Albany_ProblemUtils.cpp index 253f46dda5..0d01ae20be 100644 --- a/src/problems/Albany_ProblemUtils.cpp +++ b/src/problems/Albany_ProblemUtils.cpp @@ -134,14 +134,6 @@ getIntrepid2Basis(const CellTopologyData& ctd, bool compositeTet) return intrepidBasis; } -bool mesh_depends_on_solution () { -#ifdef ALBANY_MESH_DEPENDS_ON_SOLUTION - return true; -#else - return false; -#endif -} - bool mesh_depends_on_parameters () { #ifdef ALBANY_MESH_DEPENDS_ON_PARAMETERS return true;