MUQ  0.4.3
FlowEquation.h
Go to the documentation of this file.
1 #ifndef FLOWEQUATION_H
2 #define FLOWEQUATION_H
3 
5 
6 #include <vector>
7 
8 #include <Eigen/Core>
9 #include <Eigen/Sparse>
10 #include <Eigen/SparseCholesky>
11 
12 
13 /***
14 ## Class Definition
15 
16 This class solves the 1D elliptic PDE of the form
17  $$
18  -\frac{\partial}{\partial x}\cdot(K(x) \frac{\partial h}{\partial x}) = f(x).
19  $$
20  over $x\in[0,1]$ with boundary conditions $h(0)=0$ and
21  $\partial h/\partial x =0$ at $x=1$. This equation is a basic model of steady
22  state fluid flow in a porous media, where $h(x)$ is the hydraulic head, $K(x)$
23  is the hydraulic conductivity, and $f(x)$ is the recharge.
24 
25  This ModPiece uses linear finite elements on a uniform grid. There is a single input,
26  the conductivity $k(x)$, which is represented as piecewise constant within each
27  of the $N$ cells. There is a single output of this ModPiece: the head $h(x)$ at the
28  $N+1$ nodes in the discretization.
29 
30  INPUTS:
31 
32 
33 */
34 class FlowEquation : public muq::Modeling::ModPiece
35 {
36 public:
37 
42  FlowEquation(Eigen::VectorXd const& sourceTerm) : muq::Modeling::ModPiece({int(sourceTerm.size())},
43  {int(sourceTerm.size()+1)})
44  {
45  numCells = sourceTerm.size();
46  numNodes = sourceTerm.size()+1;
47 
48  xs = Eigen::VectorXd::LinSpaced(numNodes,0,1);
49  dx = xs(1)-xs(0);
50 
51  rhs = BuildRhs(sourceTerm);
52  };
53 
54 protected:
55 
64  void EvaluateImpl(muq::Modeling::ref_vector<Eigen::VectorXd> const& inputs) override
65  {
66  // Extract the conductivity vector from the inptus
67  auto& cond = inputs.at(0).get();
68 
69  // Build the stiffness matrix
70  auto K = BuildStiffness(cond);
71 
72  // Solve the sparse linear system and store the solution in the outputs vector
73  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
74  solver.compute(K);
75 
76  outputs.resize(1);
77  outputs.at(0) = solver.solve(rhs);
78  };
79 
108  virtual void GradientImpl(unsigned int outWrt,
109  unsigned int inWrt,
111  Eigen::VectorXd const& sens) override
112  {
113  // Extract the conductivity vector from the inptus
114  auto& cond = inputs.at(0).get();
115 
116  // Build the stiffness matrix
117  auto A = BuildStiffness(cond);
118 
119  // Factor the stiffness matrix for forward and adjoint solves
120  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
121  solver.compute(A);
122 
123  // Solve the forward problem
124  Eigen::VectorXd sol = solver.solve(rhs);
125 
126  // Solve the adjoint problem
127  Eigen::VectorXd adjRhs = BuildAdjointRHS(sens);
128  Eigen::VectorXd adjSol = solver.solve(adjRhs);
129 
130  // Compute the gradient from the adjoint solution
131  gradient = (1.0/dx)*(sol.tail(numCells)-sol.head(numCells)).array() * (adjSol.tail(numCells) - adjSol.head(numCells)).array();
132  }
133 
134 
135  /***
136  This function computes the application of the model Jacobian's matrix $J$
137  on a vector $v$. In addition to the model parameter $k$, this function also
138  requires the vector $v$.
139 
140  The gradient with respect to the conductivity field is computed by solving
141  the forward model to get $h(x)$ and then using the tangent linear approach
142  described above to obtain the Jacobian action.
143 
144  @param[in] outWrt For a model with multiple outputs, this would be the index
145  of the output list that corresponds to the sensitivity vector.
146  Since this ModPiece only has one output, the outWrt argument
147  is not used in the GradientImpl function.
148 
149  @param[in] inWrt Specifies the index of the input for which we want to compute
150  the gradient. For inWrt==0, then the Jacobian with respect
151  to the conductivity is used. Since this ModPiece only has one
152  input, 0 is the only valid value for inWrt.
153 
154  @param[in] inputs Just like the EvalauteImpl function, this is a list of vector-valued inputs.
155 
156  @param[in] vec A vector with the same size of inputs[0]. The Jacobian will be applied to this vector.
157 
158  @return This function returns nothing. It stores the result in the private
159  ModPiece::jacobianAction variable that is then returned by the `Jacobian` function.
160 
161  */
162  virtual void ApplyJacobianImpl(unsigned int outWrt,
163  unsigned int inWrt,
165  Eigen::VectorXd const& vec) override
166  {
167  // Extract the conductivity vector from the inptus
168  auto& cond = inputs.at(0).get();
169 
170  // Build the stiffness matrix
171  auto A = BuildStiffness(cond);
172 
173  // Factor the stiffness matrix for forward and tangent linear solves
174  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
175  solver.compute(A);
176 
177  // Solve the forward system
178  Eigen::VectorXd sol = solver.solve(rhs);
179 
180  // Build the tangent linear rhs
181  Eigen::VectorXd incrRhs = Eigen::VectorXd::Zero(numNodes);
182 
183  Eigen::VectorXd dh_dx = (sol.tail(numCells)-sol.head(numCells))/dx; // Spatial derivative of solution
184 
185  incrRhs.head(numCells) += ( vec.array()*dh_dx.array() ).matrix();
186  incrRhs.tail(numCells) -= ( vec.array()*dh_dx.array() ).matrix();
187 
188  // Solve the tangent linear model for the jacobian action
189  jacobianAction = solver.solve(incrRhs);
190  }
191 
230  virtual void ApplyHessianImpl(unsigned int outWrt,
231  unsigned int inWrt1,
232  unsigned int inWrt2,
234  Eigen::VectorXd const& sens,
235  Eigen::VectorXd const& vec) override
236  {
237  // Extract the conductivity vector from the inptus
238  auto& cond = inputs.at(0).get();
239 
240  // Build the stiffness matrix
241  auto K = BuildStiffness(cond);
242 
243  // Factor the stiffness matrix for forward and adjoint solves
244  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
245  solver.compute(K);
246 
247  // Solve the forward problem
248  Eigen::VectorXd sol = solver.solve(rhs);
249 
250  // If we're using the Hessian $\nabla_{kk} J$
251  if((inWrt1==0)&&(inWrt2==0)){
252 
253  // Solve the adjoint problem
254  Eigen::VectorXd adjRhs = BuildAdjointRHS(sens);
255  Eigen::VectorXd adjSol = solver.solve(adjRhs);
256 
257  // Solve the incremental forward problem
258  Eigen::VectorXd incrForRhs = BuildIncrementalRhs(sol, vec);
259  Eigen::VectorXd incrForSol = solver.solve(incrForRhs);
260 
261  // Solve the incremental adjoint problem
262  Eigen::VectorXd incrAdjRhs = BuildIncrementalRhs(adjSol, vec);
263  Eigen::VectorXd incrAdjSol = solver.solve(incrAdjRhs);
264 
265  // Construct the Hessian action
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;
270 
271  hessAction = -(incrAdjDeriv.array() * solDeriv.array() + incrForDeriv.array() * adjDeriv.array());
272 
273  // If we're using the mixed Hessian $\nabla_{ks} J$
274  }else if((inWrt1==0)&&(inWrt2==1)){
275 
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;
279 
280  hessAction = -dx * solDeriv.array() * tempDeriv.array();
281 
282  // We should never see any other options...
283  }else{
284  assert(false);
285  }
286  }
287 
289  Eigen::VectorXd BuildRhs(Eigen::VectorXd const& sourceTerm) const
290  {
291  Eigen::VectorXd rhs = Eigen::VectorXd::Zero(numNodes);
292 
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);
295 
296  return rhs;
297  }
298 
300  Eigen::VectorXd BuildAdjointRHS(Eigen::VectorXd const& sensitivity) const
301  {
302  Eigen::VectorXd rhs = -1.0*sensitivity;
303  rhs(0) = 0.0; // <- To enforce Dirichlet BC
304  return rhs;
305  }
306 
308  Eigen::VectorXd BuildIncrementalRhs(Eigen::VectorXd const& sol, Eigen::VectorXd const& khat)
309  {
310  // Compute the derivative of the solution in each cell
311  Eigen::VectorXd solGrad = (sol.tail(numCells)-sol.head(numCells))/dx;
312  Eigen::VectorXd rhs = Eigen::VectorXd::Zero(numNodes);
313 
314  unsigned int leftNode, rightNode;
315  for(unsigned int cellInd=0; cellInd<numCells; ++cellInd)
316  {
317  leftNode = cellInd;
318  rightNode = cellInd + 1;
319 
320  rhs(leftNode) -= dx*khat(cellInd)*solGrad(cellInd);
321  rhs(rightNode) += dx*khat(cellInd)*solGrad(cellInd);
322  }
323 
324  rhs(0) = 0.0; // <- To enforce Dirichlet BC at x=0
325  return rhs;
326  }
327 
329  Eigen::SparseMatrix<double> BuildStiffness(Eigen::VectorXd const& condVals) const{
330 
331  typedef Eigen::Triplet<double> T;
332  std::vector<T> nzVals;
333 
334  // Add a large number to K[0,0] to enforce Dirichlet BC
335  nzVals.push_back( T(0,0,1e10) );
336 
337  unsigned int leftNode, rightNode;
338  for(unsigned int cellInd=0; cellInd<numCells; ++cellInd)
339  {
340  leftNode = cellInd;
341  rightNode = cellInd+1;
342 
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) );
347  }
348 
349  Eigen::SparseMatrix<double> stiffMat(numNodes,numNodes);
350  stiffMat.setFromTriplets(nzVals.begin(), nzVals.end());
351 
352  return stiffMat;
353  }
354 
355 
356 private:
357 
358  // Store "mesh" information. xs contains the node locations and dx is the uniform spacing between nodes
359  Eigen::VectorXd xs;
360  double dx;
361  unsigned int numCells;
362  unsigned int numNodes;
363 
364  // Will store the precomputed RHS for the forward problem
365  Eigen::VectorXd rhs;
366 
367 }; // end of class FlowEquation
368 
369 
370 
371 
372 #endif
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)