MUQ  0.4.3
FlowEquation.h
Go to the documentation of this file.
2 
3 #include <boost/property_tree/ptree.hpp>
4 
5 #include <Eigen/Core>
6 #include <Eigen/Sparse>
7 
8 /***
9 ## Class Definition
10 
11 This class solves the 1D elliptic PDE of the form
12  $$
13  -\frac{\partial}{\partial x}\cdot(K(x) \frac{\partial h}{\partial x}) = f(x).
14  $$
15  over $x\in[0,1]$ with boundary conditions $h(0)=0$ and
16  $\partial h/\partial x =0$ at $x=1$. This equation is a basic model of steady
17  state fluid flow in a porous media, where $h(x)$ is the hydraulic head, $K(x)$
18  is the hydraulic conductivity, and $f(x)$ is the recharge.
19 
20  This ModPiece uses linear finite elements on a uniform grid. There is a single input,
21  the conductivity $k(x)$, which is represented as piecewise constant within each
22  of the $N$ cells. There is a single output of this ModPiece: the head $h(x)$ at the
23  $N+1$ nodes in the discretization.
24 
25 */
26 class FlowEquation : public muq::Modeling::ModPiece
27 {
28 public:
29 
34  FlowEquation(Eigen::VectorXd const& sourceTerm) : muq::Modeling::ModPiece({int(sourceTerm.size())},
35  {int(sourceTerm.size()+1)})
36  {
37  numCells = sourceTerm.size();
38  numNodes = sourceTerm.size()+1;
39 
40  xs = Eigen::VectorXd::LinSpaced(numNodes,0,1);
41  dx = xs(1)-xs(0);
42 
43  rhs = BuildRhs(sourceTerm);
44  };
45 
46 protected:
47 
56  void EvaluateImpl(muq::Modeling::ref_vector<Eigen::VectorXd> const& inputs) override
57  {
58  // Extract the conductivity vector from the inptus
59  auto& cond = inputs.at(0).get();
60 
61  // Build the stiffness matrix
62  auto K = BuildStiffness(cond);
63 
64  // Solve the sparse linear system and store the solution in the outputs vector
65  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
66  solver.compute(K);
67 
68  outputs.resize(1);
69  outputs.at(0) = solver.solve(rhs);
70  };
71 
100  virtual void GradientImpl(unsigned int outWrt,
101  unsigned int inWrt,
103  Eigen::VectorXd const& sens) override
104  {
105  // Extract the conductivity vector from the inptus
106  auto& cond = inputs.at(0).get();
107 
108  // Build the stiffness matrix
109  auto A = BuildStiffness(cond);
110 
111  // Factor the stiffness matrix for forward and adjoint solves
112  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
113  solver.compute(A);
114 
115  // Solve the forward problem
116  Eigen::VectorXd sol = solver.solve(rhs);
117 
118  // Solve the adjoint problem
119  Eigen::VectorXd adjRhs = BuildAdjointRHS(sens);
120  Eigen::VectorXd adjSol = solver.solve(adjRhs);
121 
122  // Compute the gradient from the adjoint solution
123  gradient = (1.0/dx)*(sol.tail(numCells)-sol.head(numCells)).array() * (adjSol.tail(numCells) - adjSol.head(numCells)).array();
124  }
125 
126 
127  /***
128  This function computes the application of the model Jacobian's matrix $J$
129  on a vector $v$. In addition to the model parameter $k$, this function also
130  requires the vector $v$.
131 
132  The gradient with respect to the conductivity field is computed by solving
133  the forward model to get $h(x)$ and then using the tangent linear approach
134  described above to obtain the Jacobian action.
135 
136  @param[in] outWrt For a model with multiple outputs, this would be the index
137  of the output list that corresponds to the sensitivity vector.
138  Since this ModPiece only has one output, the outWrt argument
139  is not used in the GradientImpl function.
140 
141  @param[in] inWrt Specifies the index of the input for which we want to compute
142  the gradient. For inWrt==0, then the Jacobian with respect
143  to the conductivity is used. Since this ModPiece only has one
144  input, 0 is the only valid value for inWrt.
145 
146  @param[in] inputs Just like the EvalauteImpl function, this is a list of vector-valued inputs.
147 
148  @param[in] vec A vector with the same size of inputs[0]. The Jacobian will be applied to this vector.
149 
150  @return This function returns nothing. It stores the result in the private
151  ModPiece::jacobianAction variable that is then returned by the `Jacobian` function.
152 
153  */
154  virtual void ApplyJacobianImpl(unsigned int outWrt,
155  unsigned int inWrt,
157  Eigen::VectorXd const& vec) override
158  {
159  // Extract the conductivity vector from the inptus
160  auto& cond = inputs.at(0).get();
161 
162  // Build the stiffness matrix
163  auto A = BuildStiffness(cond);
164 
165  // Factor the stiffness matrix for forward and tangent linear solves
166  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
167  solver.compute(A);
168 
169  // Solve the forward system
170  Eigen::VectorXd sol = solver.solve(rhs);
171 
172  // Build the tangent linear rhs
173  Eigen::VectorXd incrRhs = Eigen::VectorXd::Zero(numNodes);
174 
175  Eigen::VectorXd dh_dx = (sol.tail(numCells)-sol.head(numCells))/dx; // Spatial derivative of solution
176 
177  incrRhs.head(numCells) += ( vec.array()*dh_dx.array() ).matrix();
178  incrRhs.tail(numCells) -= ( vec.array()*dh_dx.array() ).matrix();
179 
180  // Solve the tangent linear model for the jacobian action
181  jacobianAction = solver.solve(incrRhs);
182  }
183 
222  virtual void ApplyHessianImpl(unsigned int outWrt,
223  unsigned int inWrt1,
224  unsigned int inWrt2,
226  Eigen::VectorXd const& sens,
227  Eigen::VectorXd const& vec) override
228  {
229  // Extract the conductivity vector from the inptus
230  auto& cond = inputs.at(0).get();
231 
232  // Build the stiffness matrix
233  auto K = BuildStiffness(cond);
234 
235  // Factor the stiffness matrix for forward and adjoint solves
236  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
237  solver.compute(K);
238 
239  // Solve the forward problem
240  Eigen::VectorXd sol = solver.solve(rhs);
241 
242  // If we're using the Hessian $\nabla_{kk} J$
243  if((inWrt1==0)&&(inWrt2==0)){
244 
245  // Solve the adjoint problem
246  Eigen::VectorXd adjRhs = BuildAdjointRHS(sens);
247  Eigen::VectorXd adjSol = solver.solve(adjRhs);
248 
249  // Solve the incremental forward problem
250  Eigen::VectorXd incrForRhs = BuildIncrementalRhs(sol, vec);
251  Eigen::VectorXd incrForSol = solver.solve(incrForRhs);
252 
253  // Solve the incremental adjoint problem
254  Eigen::VectorXd incrAdjRhs = BuildIncrementalRhs(adjSol, vec);
255  Eigen::VectorXd incrAdjSol = solver.solve(incrAdjRhs);
256 
257  // Construct the Hessian action
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;
262 
263  hessAction = -(incrAdjDeriv.array() * solDeriv.array() + incrForDeriv.array() * adjDeriv.array());
264 
265  // If we're using the mixed Hessian $\nabla_{ks} J$
266  }else if((inWrt1==0)&&(inWrt2==1)){
267 
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;
271 
272  hessAction = -dx * solDeriv.array() * tempDeriv.array();
273 
274  // We should never see any other options...
275  }else{
276  assert(false);
277  }
278  }
279 
281  Eigen::VectorXd BuildRhs(Eigen::VectorXd const& sourceTerm) const
282  {
283  Eigen::VectorXd rhs = Eigen::VectorXd::Zero(numNodes);
284 
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);
287 
288  return rhs;
289  }
290 
292  Eigen::VectorXd BuildAdjointRHS(Eigen::VectorXd const& sensitivity) const
293  {
294  Eigen::VectorXd rhs = -1.0*sensitivity;
295  rhs(0) = 0.0; // <- To enforce Dirichlet BC
296  return rhs;
297  }
298 
300  Eigen::VectorXd BuildIncrementalRhs(Eigen::VectorXd const& sol, Eigen::VectorXd const& khat)
301  {
302  // Compute the derivative of the solution in each cell
303  Eigen::VectorXd solGrad = (sol.tail(numCells)-sol.head(numCells))/dx;
304  Eigen::VectorXd rhs = Eigen::VectorXd::Zero(numNodes);
305 
306  unsigned int leftNode, rightNode;
307  for(unsigned int cellInd=0; cellInd<numCells; ++cellInd)
308  {
309  leftNode = cellInd;
310  rightNode = cellInd + 1;
311 
312  rhs(leftNode) -= dx*khat(cellInd)*solGrad(cellInd);
313  rhs(rightNode) += dx*khat(cellInd)*solGrad(cellInd);
314  }
315 
316  rhs(0) = 0.0; // <- To enforce Dirichlet BC at x=0
317  return rhs;
318  }
319 
321  Eigen::SparseMatrix<double> BuildStiffness(Eigen::VectorXd const& condVals) const{
322 
323  typedef Eigen::Triplet<double> T;
324  std::vector<T> nzVals;
325 
326  // Add a large number to K[0,0] to enforce Dirichlet BC
327  nzVals.push_back( T(0,0,1e10) );
328 
329  unsigned int leftNode, rightNode;
330  for(unsigned int cellInd=0; cellInd<numCells; ++cellInd)
331  {
332  leftNode = cellInd;
333  rightNode = cellInd+1;
334 
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) );
339  }
340 
341  Eigen::SparseMatrix<double> stiffMat(numNodes,numNodes);
342  stiffMat.setFromTriplets(nzVals.begin(), nzVals.end());
343 
344  return stiffMat;
345  }
346 
347 
348 private:
349 
350  // Store "mesh" information. xs contains the node locations and dx is the uniform spacing between nodes
351  Eigen::VectorXd xs;
352  double dx;
353  unsigned int numCells;
354  unsigned int numNodes;
355 
356  // Will store the precomputed RHS for the forward problem
357  Eigen::VectorXd rhs;
358 
359 }; // end of class SimpleModel
Provides an abstract interface for defining vector-valued model components.
Definition: ModPiece.h:148
virtual void GradientImpl(unsigned int const outputDimWrt, unsigned int const inputDimWrt, ref_vector< Eigen::VectorXd > const &input, Eigen::VectorXd const &sensitivity)
Definition: ModPiece.cpp:237
virtual void ApplyJacobianImpl(unsigned int const outputDimWrt, unsigned int const inputDimWrt, ref_vector< Eigen::VectorXd > const &input, Eigen::VectorXd const &vec)
Definition: ModPiece.cpp:262
std::vector< Eigen::VectorXd > outputs
Definition: ModPiece.h:503
ModPiece(std::vector< int > const &inputSizes, std::vector< int > const &outputSizes)
Definition: ModPiece.cpp:9
virtual void EvaluateImpl(ref_vector< boost::any > const &inputs) override
User-implemented function that determines the behavior of this muq::Modeling::WorkPiece.
Definition: ModPiece.cpp:179
Eigen::VectorXd gradient
Definition: ModPiece.h:504
Eigen::VectorXd hessAction
Definition: ModPiece.h:507
Eigen::VectorXd jacobianAction
Definition: ModPiece.h:505
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)
Definition: ModPiece.cpp:436
std::vector< std::reference_wrapper< const T > > ref_vector
A vector of references to something ...
Definition: WorkPiece.h:37
@ array
array (ordered collection of values)