3 #include <boost/property_tree/ptree.hpp>
6 #include <Eigen/Sparse>
34 FlowEquation(Eigen::VectorXd
const& sourceTerm) :
muq::Modeling::
ModPiece({int(sourceTerm.size())},
35 {int(sourceTerm.size()+1)})
37 numCells = sourceTerm.size();
38 numNodes = sourceTerm.size()+1;
40 xs = Eigen::VectorXd::LinSpaced(numNodes,0,1);
43 rhs = BuildRhs(sourceTerm);
59 auto& cond = inputs.at(0).get();
62 auto K = BuildStiffness(cond);
65 Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
69 outputs.at(0) = solver.solve(rhs);
103 Eigen::VectorXd
const& sens)
override
106 auto& cond = inputs.at(0).get();
109 auto A = BuildStiffness(cond);
112 Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
116 Eigen::VectorXd sol = solver.solve(rhs);
119 Eigen::VectorXd adjRhs = BuildAdjointRHS(sens);
120 Eigen::VectorXd adjSol = solver.solve(adjRhs);
123 gradient = (1.0/dx)*(sol.tail(numCells)-sol.head(numCells)).
array() * (adjSol.tail(numCells) - adjSol.head(numCells)).
array();
157 Eigen::VectorXd
const& vec)
override
160 auto& cond = inputs.at(0).get();
163 auto A = BuildStiffness(cond);
166 Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
170 Eigen::VectorXd sol = solver.solve(rhs);
173 Eigen::VectorXd incrRhs = Eigen::VectorXd::Zero(numNodes);
175 Eigen::VectorXd dh_dx = (sol.tail(numCells)-sol.head(numCells))/dx;
177 incrRhs.head(numCells) += ( vec.array()*dh_dx.array() ).matrix();
178 incrRhs.tail(numCells) -= ( vec.array()*dh_dx.array() ).matrix();
226 Eigen::VectorXd
const& sens,
227 Eigen::VectorXd
const& vec)
override
230 auto& cond = inputs.at(0).get();
233 auto K = BuildStiffness(cond);
236 Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
240 Eigen::VectorXd sol = solver.solve(rhs);
243 if((inWrt1==0)&&(inWrt2==0)){
246 Eigen::VectorXd adjRhs = BuildAdjointRHS(sens);
247 Eigen::VectorXd adjSol = solver.solve(adjRhs);
250 Eigen::VectorXd incrForRhs = BuildIncrementalRhs(sol, vec);
251 Eigen::VectorXd incrForSol = solver.solve(incrForRhs);
254 Eigen::VectorXd incrAdjRhs = BuildIncrementalRhs(adjSol, vec);
255 Eigen::VectorXd incrAdjSol = solver.solve(incrAdjRhs);
258 auto solDeriv = (sol.tail(numCells)-sol.head(numCells))/dx;
259 auto adjDeriv = (adjSol.tail(numCells)-adjSol.head(numCells))/dx;
260 auto incrForDeriv = (incrForSol.tail(numCells) - incrForSol.head(numCells))/dx;
261 auto incrAdjDeriv = (incrAdjSol.tail(numCells) - incrAdjSol.head(numCells))/dx;
263 hessAction = -(incrAdjDeriv.array() * solDeriv.array() + incrForDeriv.array() * adjDeriv.array());
266 }
else if((inWrt1==0)&&(inWrt2==1)){
268 Eigen::VectorXd temp = solver.solve(vec);
269 auto solDeriv = (sol.tail(numCells) - sol.head(numCells))/dx;
270 auto tempDeriv = (temp.tail(numCells)-temp.head(numCells))/dx;
272 hessAction = -dx * solDeriv.array() * tempDeriv.array();
281 Eigen::VectorXd BuildRhs(Eigen::VectorXd
const& sourceTerm)
const
283 Eigen::VectorXd rhs = Eigen::VectorXd::Zero(numNodes);
285 rhs.segment(1,numNodes-2) = 0.5*dx*(sourceTerm.tail(numNodes-2) + sourceTerm.head(numNodes-2));
286 rhs(numNodes-1) = 0.5*dx*sourceTerm(numNodes-2);
292 Eigen::VectorXd BuildAdjointRHS(Eigen::VectorXd
const& sensitivity)
const
294 Eigen::VectorXd rhs = -1.0*sensitivity;
300 Eigen::VectorXd BuildIncrementalRhs(Eigen::VectorXd
const& sol, Eigen::VectorXd
const& khat)
303 Eigen::VectorXd solGrad = (sol.tail(numCells)-sol.head(numCells))/dx;
304 Eigen::VectorXd rhs = Eigen::VectorXd::Zero(numNodes);
306 unsigned int leftNode, rightNode;
307 for(
unsigned int cellInd=0; cellInd<numCells; ++cellInd)
310 rightNode = cellInd + 1;
312 rhs(leftNode) -= dx*khat(cellInd)*solGrad(cellInd);
313 rhs(rightNode) += dx*khat(cellInd)*solGrad(cellInd);
321 Eigen::SparseMatrix<double> BuildStiffness(Eigen::VectorXd
const& condVals)
const{
323 typedef Eigen::Triplet<double> T;
324 std::vector<T> nzVals;
327 nzVals.push_back( T(0,0,1e10) );
329 unsigned int leftNode, rightNode;
330 for(
unsigned int cellInd=0; cellInd<numCells; ++cellInd)
333 rightNode = cellInd+1;
335 nzVals.push_back( T(leftNode, rightNode, -condVals(cellInd)/dx) );
336 nzVals.push_back( T(rightNode, leftNode, -condVals(cellInd)/dx) );
337 nzVals.push_back( T(rightNode, rightNode, condVals(cellInd)/dx) );
338 nzVals.push_back( T(leftNode, leftNode, condVals(cellInd)/dx) );
341 Eigen::SparseMatrix<double> stiffMat(numNodes,numNodes);
342 stiffMat.setFromTriplets(nzVals.begin(), nzVals.end());
353 unsigned int numCells;
354 unsigned int numNodes;
Provides an abstract interface for defining vector-valued model components.
virtual void GradientImpl(unsigned int const outputDimWrt, unsigned int const inputDimWrt, ref_vector< Eigen::VectorXd > const &input, Eigen::VectorXd const &sensitivity)
virtual void ApplyJacobianImpl(unsigned int const outputDimWrt, unsigned int const inputDimWrt, ref_vector< Eigen::VectorXd > const &input, Eigen::VectorXd const &vec)
std::vector< Eigen::VectorXd > outputs
ModPiece(std::vector< int > const &inputSizes, std::vector< int > const &outputSizes)
virtual void EvaluateImpl(ref_vector< boost::any > const &inputs) override
User-implemented function that determines the behavior of this muq::Modeling::WorkPiece.
Eigen::VectorXd hessAction
Eigen::VectorXd jacobianAction
virtual void ApplyHessianImpl(unsigned int const outWrt, unsigned int const inWrt1, unsigned int const inWrt2, ref_vector< Eigen::VectorXd > const &input, Eigen::VectorXd const &sens, Eigen::VectorXd const &vec)
std::vector< std::reference_wrapper< const T > > ref_vector
A vector of references to something ...
@ array
array (ordered collection of values)