9 #include <Eigen/Sparse>
10 #include <Eigen/SparseCholesky>
42 FlowEquation(Eigen::VectorXd
const& sourceTerm) :
muq::Modeling::
ModPiece({int(sourceTerm.size())},
43 {int(sourceTerm.size()+1)})
45 numCells = sourceTerm.size();
46 numNodes = sourceTerm.size()+1;
48 xs = Eigen::VectorXd::LinSpaced(numNodes,0,1);
51 rhs = BuildRhs(sourceTerm);
67 auto& cond = inputs.at(0).get();
70 auto K = BuildStiffness(cond);
73 Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
77 outputs.at(0) = solver.solve(rhs);
111 Eigen::VectorXd
const& sens)
override
114 auto& cond = inputs.at(0).get();
117 auto A = BuildStiffness(cond);
120 Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
124 Eigen::VectorXd sol = solver.solve(rhs);
127 Eigen::VectorXd adjRhs = BuildAdjointRHS(sens);
128 Eigen::VectorXd adjSol = solver.solve(adjRhs);
131 gradient = (1.0/dx)*(sol.tail(numCells)-sol.head(numCells)).
array() * (adjSol.tail(numCells) - adjSol.head(numCells)).
array();
165 Eigen::VectorXd
const& vec)
override
168 auto& cond = inputs.at(0).get();
171 auto A = BuildStiffness(cond);
174 Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
178 Eigen::VectorXd sol = solver.solve(rhs);
181 Eigen::VectorXd incrRhs = Eigen::VectorXd::Zero(numNodes);
183 Eigen::VectorXd dh_dx = (sol.tail(numCells)-sol.head(numCells))/dx;
185 incrRhs.head(numCells) += ( vec.array()*dh_dx.array() ).matrix();
186 incrRhs.tail(numCells) -= ( vec.array()*dh_dx.array() ).matrix();
234 Eigen::VectorXd
const& sens,
235 Eigen::VectorXd
const& vec)
override
238 auto& cond = inputs.at(0).get();
241 auto K = BuildStiffness(cond);
244 Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
248 Eigen::VectorXd sol = solver.solve(rhs);
251 if((inWrt1==0)&&(inWrt2==0)){
254 Eigen::VectorXd adjRhs = BuildAdjointRHS(sens);
255 Eigen::VectorXd adjSol = solver.solve(adjRhs);
258 Eigen::VectorXd incrForRhs = BuildIncrementalRhs(sol, vec);
259 Eigen::VectorXd incrForSol = solver.solve(incrForRhs);
262 Eigen::VectorXd incrAdjRhs = BuildIncrementalRhs(adjSol, vec);
263 Eigen::VectorXd incrAdjSol = solver.solve(incrAdjRhs);
266 auto solDeriv = (sol.tail(numCells)-sol.head(numCells))/dx;
267 auto adjDeriv = (adjSol.tail(numCells)-adjSol.head(numCells))/dx;
268 auto incrForDeriv = (incrForSol.tail(numCells) - incrForSol.head(numCells))/dx;
269 auto incrAdjDeriv = (incrAdjSol.tail(numCells) - incrAdjSol.head(numCells))/dx;
271 hessAction = -(incrAdjDeriv.array() * solDeriv.array() + incrForDeriv.array() * adjDeriv.array());
274 }
else if((inWrt1==0)&&(inWrt2==1)){
276 Eigen::VectorXd temp = solver.solve(vec);
277 auto solDeriv = (sol.tail(numCells) - sol.head(numCells))/dx;
278 auto tempDeriv = (temp.tail(numCells)-temp.head(numCells))/dx;
280 hessAction = -dx * solDeriv.array() * tempDeriv.array();
289 Eigen::VectorXd BuildRhs(Eigen::VectorXd
const& sourceTerm)
const
291 Eigen::VectorXd rhs = Eigen::VectorXd::Zero(numNodes);
293 rhs.segment(1,numNodes-2) = 0.5*dx*(sourceTerm.tail(numNodes-2) + sourceTerm.head(numNodes-2));
294 rhs(numNodes-1) = 0.5*dx*sourceTerm(numNodes-2);
300 Eigen::VectorXd BuildAdjointRHS(Eigen::VectorXd
const& sensitivity)
const
302 Eigen::VectorXd rhs = -1.0*sensitivity;
308 Eigen::VectorXd BuildIncrementalRhs(Eigen::VectorXd
const& sol, Eigen::VectorXd
const& khat)
311 Eigen::VectorXd solGrad = (sol.tail(numCells)-sol.head(numCells))/dx;
312 Eigen::VectorXd rhs = Eigen::VectorXd::Zero(numNodes);
314 unsigned int leftNode, rightNode;
315 for(
unsigned int cellInd=0; cellInd<numCells; ++cellInd)
318 rightNode = cellInd + 1;
320 rhs(leftNode) -= dx*khat(cellInd)*solGrad(cellInd);
321 rhs(rightNode) += dx*khat(cellInd)*solGrad(cellInd);
329 Eigen::SparseMatrix<double> BuildStiffness(Eigen::VectorXd
const& condVals)
const{
331 typedef Eigen::Triplet<double> T;
332 std::vector<T> nzVals;
335 nzVals.push_back( T(0,0,1e10) );
337 unsigned int leftNode, rightNode;
338 for(
unsigned int cellInd=0; cellInd<numCells; ++cellInd)
341 rightNode = cellInd+1;
343 nzVals.push_back( T(leftNode, rightNode, -condVals(cellInd)/dx) );
344 nzVals.push_back( T(rightNode, leftNode, -condVals(cellInd)/dx) );
345 nzVals.push_back( T(rightNode, rightNode, condVals(cellInd)/dx) );
346 nzVals.push_back( T(leftNode, leftNode, condVals(cellInd)/dx) );
349 Eigen::SparseMatrix<double> stiffMat(numNodes,numNodes);
350 stiffMat.setFromTriplets(nzVals.begin(), nzVals.end());
361 unsigned int numCells;
362 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)