Skip to content

Commit

Permalink
Changing how parameter for sensitivity analysis (kappa_x for now)
Browse files Browse the repository at this point in the history 
is registered.

Registering it as a shared parameter (LandIce::SharedParameter) rather than
Sacado parameter.  I'm not sure why the Sacado parameter doesn't work, but I
suppose it doesn't really matter since we will be using LandIce::SharedParameter
for all intents and purposes.

With this change, dg/dp \neq 0 for the Solution Average response, and dg/dp = 0
for the Solution Max Value response, as expected.  I have not verified correctness
for the former response yet.  Before, I was always getting df/dp = 0 and dg/dp = 0
when attempting to compute transient sensitivities w.r.t. kappa_x.

Code needs some cleanup - exception checks, ifdefs, and extension to more
parameters (kappa_y and kappa_z).
  • Loading branch information
@ikalash
ikalash committed May 27, 2020
1 parent 5fd814a commit a127d01
Show file tree
Hide file tree
Showing 4 changed files with 49 additions and 63 deletions.
13 changes: 7 additions & 6 deletions src/evaluators/pde/PHAL_ThermalResidWithSensitivities.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include "Phalanx_Evaluator_Derived.hpp"
#include "Phalanx_MDField.hpp"
#include "Sacado_ParameterAccessor.hpp"
#include "Albany_Layouts.hpp"

namespace PHAL {

Expand All @@ -23,12 +24,12 @@ namespace PHAL {

template<typename EvalT, typename Traits>
class ThermalResidWithSensitivities : public PHX::EvaluatorWithBaseImpl<Traits>,
public PHX::EvaluatorDerived<EvalT, Traits>,
public Sacado::ParameterAccessor<EvalT, SPL_Traits> {
public PHX::EvaluatorDerived<EvalT, Traits> {

public:

ThermalResidWithSensitivities(Teuchos::ParameterList const& p);
ThermalResidWithSensitivities(Teuchos::ParameterList const& p,
const Teuchos::RCP<Albany::Layouts>& dl);

void
postRegistrationSetup(
Expand All @@ -38,8 +39,6 @@ class ThermalResidWithSensitivities : public PHX::EvaluatorWithBaseImpl<Traits>,
void
evaluateFields(typename Traits::EvalData d);

typename EvalT::ScalarT& getValue(const std::string &n);

private:

typedef typename EvalT::ScalarT ScalarT;
Expand All @@ -53,9 +52,11 @@ class ThermalResidWithSensitivities : public PHX::EvaluatorWithBaseImpl<Traits>,
PHX::MDField<const MeshScalarT, Cell, Node, QuadPoint, Dim> wGradBF;
PHX::MDField<ScalarT const, Cell, QuadPoint, Dim> TGrad;
PHX::MDField<const MeshScalarT, Cell, QuadPoint, Dim> coordVec;
Teuchos::Array<double> kappa; // Thermal Conductivity array
double C; // Heat Capacity
double rho; // Density
PHX::MDField<const ScalarT> kappa_x;
double kappa_y;
double kappa_z;

// Output:
PHX::MDField<ScalarT, Cell, QuadPoint> Source;
Expand Down
73 changes: 23 additions & 50 deletions src/evaluators/pde/PHAL_ThermalResidWithSensitivities_Def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,8 @@ namespace PHAL {
//**********************************************************************
template<typename EvalT, typename Traits>
ThermalResidWithSensitivities<EvalT, Traits>::
ThermalResidWithSensitivities(const Teuchos::ParameterList& p) :
ThermalResidWithSensitivities(const Teuchos::ParameterList& p,
const Teuchos::RCP<Albany::Layouts>& dl) :
wBF (p.get<std::string> ("Weighted BF Name"),
p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Scalar Data Layout") ),
Tdot (p.get<std::string> ("QP Time Derivative Variable Name"),
Expand All @@ -30,15 +31,18 @@ ThermalResidWithSensitivities(const Teuchos::ParameterList& p) :
p.get<Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout") ),
Source (p.get<std::string> ("Source Name"),
p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
kappa(p.get<Teuchos::Array<double>>("Thermal Conductivity")),
rho(p.get<double>("Density")),
C(p.get<double>("Heat Capacity"))
C(p.get<double>("Heat Capacity")),
kappa_x("kappa_x Parameter", dl->shared_param),
kappa_y(p.get<double>("kappa_y")),
kappa_z(p.get<double>("kappa_z"))
{

this->addDependentField(wBF);
this->addDependentField(Tdot);
this->addDependentField(TGrad);
this->addDependentField(wGradBF);
this->addDependentField(kappa_x);
this->addEvaluatedField(Source);
this->addEvaluatedField(TResidual);
Teuchos::RCP<PHX::DataLayout> vector_dl =
Expand All @@ -55,12 +59,13 @@ ThermalResidWithSensitivities(const Teuchos::ParameterList& p) :
numQPs = dims[2];
numDims = dims[3];

if (kappa.size() != numDims) {
//IKT FIXME: add error checking to ensure size(kappa) == numDims
/*if (kappa.size() != numDims) {
TEUCHOS_TEST_FOR_EXCEPTION(
true,
std::logic_error,
"Thermal Conductivity size " << kappa.size() << " != # dimensions " << numDims << "\n");
}
}*/

std::string thermal_source = p.get<std::string>("Thermal Source");
if (thermal_source == "None") {
Expand Down Expand Up @@ -93,20 +98,6 @@ ThermalResidWithSensitivities(const Teuchos::ParameterList& p) :
"Unknown Thermal Source = " << thermal_source << "! Valid options are: 'None', '1D Cost' and '2D Cost Expt'. \n");
}

// Add kappa[0] wavelength as a Sacado-ized parameter
Teuchos::RCP<ParamLib> paramLib =
p.get< Teuchos::RCP<ParamLib> >("Parameter Library");
this->registerSacadoParameter("kappa_x", paramLib);
if (numDims > 1) {
// Add kappa[1] wavelength as a Sacado-ized parameter
this->registerSacadoParameter("kappa_y", paramLib);
}
if (numDims > 2) {
// Add kappa[2] wavelength as a Sacado-ized parameter
this->registerSacadoParameter("kappa_z", paramLib);
}
this->registerSacadoParameter("C", paramLib);
this->registerSacadoParameter("rho", paramLib);
this->setName("ThermalResidWithSensitivities" );
}

Expand All @@ -123,28 +114,7 @@ postRegistrationSetup(typename Traits::SetupData d,
this->utils.setFieldData(Source, fm);
this->utils.setFieldData(coordVec, fm);
this->utils.setFieldData(TResidual, fm);
}
// **********************************************************************
template<typename EvalT,typename Traits>
typename ThermalResidWithSensitivities<EvalT,Traits>::ScalarT&
ThermalResidWithSensitivities<EvalT,Traits>::getValue(const std::string &n)
{
if (n == "kappa_x") {
value = kappa[0];
}
else if (n == "kappa_y") {
value = kappa[1];
}
else if (n == "kappa_z") {
value = kappa[2];
}
else if (n == "C") {
value = C;
}
else if (n == "rho") {
value = rho;
}
return value;
this->utils.setFieldData(kappa_x, fm);
}

//**********************************************************************
Expand All @@ -160,27 +130,27 @@ evaluateFields(typename Traits::EvalData workset)
}
}
}
else if (force_type == ONEDCOST) { //Source term such that T = a*x*(1-x)*cos(2*pi*kappa[0]*t/rho/C) with a = 16.0
else if (force_type == ONEDCOST) { //Source term such that T = a*x*(1-x)*cos(2*pi*kappa_x*t/rho/C) with a = 16.0
for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
for (std::size_t qp=0; qp < numQPs; ++qp) {
const ScalarT x = coordVec(cell, qp, 0);
const RealType a = 16.0;
const RealType t = workset.current_time;
Source(cell, qp) = -2.0*a*kappa[0]*(M_PI*x*(1-x)*sin(2.0*M_PI*kappa[0]*t/rho/C) - cos(2.0*M_PI*kappa[0]*t/rho/C));
Source(cell, qp) = -2.0*a*kappa_x(0)*(M_PI*x*(1-x)*sin(2.0*M_PI*kappa_x(0)*t/rho/C) - cos(2.0*M_PI*kappa_x(0)*t/rho/C));
}
}
}
else if (force_type == TWODCOSTEXPT) { //Source term such that T = a*x*(1-x)*y*(1-y)*
//cos(2*pi*kappa[0]*t/rho/C)*exp(2*pi*kappa[1]*t/rho/C) with a = 16.0
//cos(2*pi*kappa_x*t/rho/C)*exp(2*pi*kappa_y*t/rho/C) with a = 16.0
for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
for (std::size_t qp=0; qp < numQPs; ++qp) {
const ScalarT x = coordVec(cell, qp, 0);
const ScalarT y = coordVec(cell, qp, 1);
const RealType a = 16.0;
const RealType t = workset.current_time;
Source(cell, qp) = 2.0*M_PI*a*x*(1.0-x)*y*(1.0-y)*exp(2.0*M_PI*kappa[1]*t/rho/C)*(kappa[1]*cos(2*M_PI*kappa[0]*t/rho/C)
-kappa[0]*sin(2.0*M_PI*kappa[0]*t/rho/C)) + 2.0*a*cos(2.0*M_PI*kappa[0]*t/rho/C)
*exp(2.0*M_PI*kappa[1]*t/rho/C)*(kappa[0]*y*(1-y) + kappa[1]*x*(1-x));
Source(cell, qp) = 2.0*M_PI*a*x*(1.0-x)*y*(1.0-y)*exp(2.0*M_PI*kappa_y*t/rho/C)*(kappa_y*cos(2*M_PI*kappa_x(0)*t/rho/C)
-kappa_x(0)*sin(2.0*M_PI*kappa_x(0)*t/rho/C)) + 2.0*a*cos(2.0*M_PI*kappa_x(0)*t/rho/C)
*exp(2.0*M_PI*kappa_y*t/rho/C)*(kappa_x(0)*y*(1-y) + kappa_y*x*(1-x));
}
}
}
Expand All @@ -197,9 +167,12 @@ evaluateFields(typename Traits::EvalData workset)
// Source contribution to residual
- Source(cell,qp) * wBF(cell, node, qp);
// Diffusion part of residual
for (std::size_t ndim = 0; ndim < numDims; ++ndim) {
TResidual(cell, node) += kappa[ndim] * TGrad(cell, qp, ndim) *
wGradBF(cell, node, qp, ndim);
TResidual(cell, node) += kappa_x(0) * TGrad(cell, qp, 0) * wGradBF(cell, node, qp, 0);
if (numDims > 1) {
TResidual(cell, node) += kappa_y * TGrad(cell, qp, 1) * wGradBF(cell, node, qp, 1);
}
if (numDims > 2) {
TResidual(cell, node) += kappa_z * TGrad(cell, qp, 2) * wGradBF(cell, node, qp, 2);
}
}
}
Expand Down
24 changes: 18 additions & 6 deletions src/problems/Albany_ThermalProblemWithSensitivities.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@

#include "PHAL_ConvertFieldType.hpp"
#include "Albany_MaterialDatabase.hpp"

namespace Albany {

/*!
Expand Down Expand Up @@ -115,7 +114,9 @@ namespace Albany {
#include "Albany_ResponseUtilities.hpp"
//#include "PHAL_Neumann.hpp"
#include "PHAL_ThermalResidWithSensitivities.hpp"

//IKT FIXME - add dependence on LandIce
#include "../../LandIce/evaluators/LandIce_SharedParameter.hpp"
#include "../../LandIce/problems/LandIce_ParamEnum.hpp"

template <typename EvalT>
Teuchos::RCP<const PHX::FieldTag>
Expand Down Expand Up @@ -204,6 +205,17 @@ Albany::ThermalProblemWithSensitivities::constructEvaluators(
evalUtils.constructDOFGradInterpolationEvaluator(dof_names[i]));
}

{ //Shared parameters for sensitivity analysis
RCP<ParameterList> p = rcp(new ParameterList("Thermal Conductivity Parameters"));
p->set< RCP<ParamLib> >("Parameter Library", paramLib);
const std::string param_name = "kappa_x Parameter";
p->set<std::string>("Parameter Name", param_name);
RCP<LandIce::SharedParameter<EvalT,PHAL::AlbanyTraits,LandIce::ParamEnum,LandIce::ParamEnum::Homotopy>> ptr_kappa_x;
ptr_kappa_x = rcp(new LandIce::SharedParameter<EvalT,PHAL::AlbanyTraits,LandIce::ParamEnum,LandIce::ParamEnum::Homotopy>(*p,dl));
ptr_kappa_x->setNominalValue(params->sublist("Parameters"), kappa[0]);
fm0.template registerEvaluator<EvalT>(ptr_kappa_x);
}

{ // Temperature Resid
RCP<ParameterList> p = rcp(new ParameterList("Temperature Resid"));

Expand All @@ -221,17 +233,17 @@ Albany::ThermalProblemWithSensitivities::constructEvaluators(
p->set<RCP<DataLayout>>("QP Vector Data Layout", dl->qp_vector);
p->set<RCP<DataLayout>>("Node QP Vector Data Layout", dl->node_qp_vector);
p->set<RCP<DataLayout>>("Node Scalar Data Layout", dl->node_scalar);

p->set<RCP<ParamLib> >("Parameter Library", paramLib);
p->set<Teuchos::Array<double>>("Thermal Conductivity", kappa);
p->set<std::string>("Continuation Parameter Name","kappa_x Parameter");
p->set<double>("Heat Capacity", C);
p->set<double>("Density", rho);
p->set<double>("kappa_y", kappa[1]);
p->set<double>("kappa_z", kappa[2]);
p->set<std::string>("Thermal Source", thermal_source);

// Output
p->set<string>("Residual Name", "Temperature Residual");

ev = rcp(new PHAL::ThermalResidWithSensitivities<EvalT, AlbanyTraits>(*p));
ev = rcp(new PHAL::ThermalResidWithSensitivities<EvalT, AlbanyTraits>(*p, dl));
fm0.template registerEvaluator<EvalT>(ev);
}

Expand Down
2 changes: 1 addition & 1 deletion tests/small/Thermal1D/input_with_source.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ ALBANY:
Response 0: Solution Average
Parameters:
Number: 1
Parameter 0: 'kappa_x'
Parameter 0: 'kappa_x Parameter'
Discretization:
1D Elements: 20
1D Scale: 1.00000000000000000e+00
Expand Down

0 comments on commit a127d01

Please sign in to comment.