Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
twoD_diffusion_problem.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
43
45#include "EpetraExt_MatrixMatrix.h"
46#include "Teuchos_Assert.hpp"
48
49namespace {
50
51 // Functor representing a diffusion function given by a KL expansion
52 struct KL_Diffusion_Func {
53 double mean;
54 mutable Teuchos::Array<double> point;
55 Teuchos::RCP< Stokhos::KL::ExponentialRandomField<double> > rf;
56
57 KL_Diffusion_Func(double xyLeft, double xyRight,
58 double mean_, double std_dev,
59 double L, int num_KL) : mean(mean_), point(2)
60 {
61 Teuchos::ParameterList rfParams;
62 rfParams.set("Number of KL Terms", num_KL);
63 rfParams.set("Mean", mean);
64 rfParams.set("Standard Deviation", std_dev);
65 int ndim = 2;
66 Teuchos::Array<double> domain_upper(ndim), domain_lower(ndim),
67 correlation_length(ndim);
68 for (int i=0; i<ndim; i++) {
69 domain_upper[i] = xyRight;
70 domain_lower[i] = xyLeft;
71 correlation_length[i] = L;
72 }
73 rfParams.set("Domain Upper Bounds", domain_upper);
74 rfParams.set("Domain Lower Bounds", domain_lower);
75 rfParams.set("Correlation Lengths", correlation_length);
76
77 rf =
78 Teuchos::rcp(new Stokhos::KL::ExponentialRandomField<double>(rfParams));
79 }
80
81 double operator() (double x, double y, int k) const {
82 if (k == 0)
83 return mean;
84 point[0] = x;
85 point[1] = y;
86 return rf->evaluate_eigenfunction(point, k-1);
87 }
88 };
89
90 // Functor representing a diffusion function given by a KL expansion
91 // where the basis is normalized
92 template <typename DiffusionFunc>
93 struct Normalized_KL_Diffusion_Func {
94 const DiffusionFunc& func;
95 int d;
96 Teuchos::Array<double> psi_0, psi_1;
97
98 Normalized_KL_Diffusion_Func(
99 const DiffusionFunc& func_,
101 func(func_),
102 d(basis.dimension()),
103 psi_0(basis.size()),
104 psi_1(basis.size())
105 {
106 Teuchos::Array<double> zero(d), one(d);
107 for(int k=0; k<d; k++) {
108 zero[k] = 0.0;
109 one[k] = 1.0;
110 }
111 basis.evaluateBases(zero, psi_0);
112 basis.evaluateBases(one, psi_1);
113 }
114
115 double operator() (double x, double y, int k) const {
116 if (k == 0) {
117 double val = func(x, y, 0);
118 for (int i=1; i<=d; i++)
119 val -= psi_0[i]/(psi_1[i]-psi_0[i])*func(x, y, i);
120 val /= psi_0[0];
121 return val;
122 }
123 else
124 return 1.0/(psi_1[k]-psi_0[k])*func(x, y, k);
125 }
126 };
127
128 // Functor representing a diffusion function given by a log-normal PC expansion
129 template <typename DiffusionFunc>
130 struct LogNormal_Diffusion_Func {
131 double mean;
132 const DiffusionFunc& func;
133 Teuchos::RCP<const Stokhos::ProductBasis<int, double> > prodbasis;
134
135 LogNormal_Diffusion_Func(
136 double mean_,
137 const DiffusionFunc& func_,
138 const Teuchos::RCP<const Stokhos::ProductBasis<int, double> > prodbasis_)
139 : mean(mean_), func(func_), prodbasis(prodbasis_) {}
140
141 double operator() (double x, double y, int k) const {
142 int d = prodbasis->dimension();
143 const Teuchos::Array<double>& norms = prodbasis->norm_squared();
144 Teuchos::Array<int> multiIndex = prodbasis->getTerm(k);
145 double sum_g = 0.0, efval;
146 for (int l=0; l<d; l++) {
147 sum_g += std::pow(func(x,y,l+1),2);
148 }
149 efval = std::exp(mean + 0.5*sum_g)/norms[k];
150 for (int l=0; l<d; l++) {
151 efval *= std::pow(func(x,y,l+1),multiIndex[l]);
152 }
153 return efval;
154 }
155 };
156
157}
158
161 const Teuchos::RCP<const Epetra_Comm>& comm, int n, int d,
162 double s, double mu,
163 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& basis_,
164 bool log_normal_,
165 bool eliminate_bcs_) :
166 mesh(n*n),
167 basis(basis_),
168 log_normal(log_normal_),
169 eliminate_bcs(eliminate_bcs_)
170{
172 // Construct the mesh.
173 // The mesh is uniform and the nodes are numbered
174 // LEFT to RIGHT, DOWN to UP.
175 //
176 // 5-6-7-8-9
177 // | | | | |
178 // 0-1-2-3-4
180 double xyLeft = -.5;
181 double xyRight = .5;
182 h = (xyRight - xyLeft)/((double)(n-1));
183 Teuchos::Array<int> global_dof_indices;
184 for (int j=0; j<n; j++) {
185 double y = xyLeft + j*h;
186 for (int i=0; i<n; i++) {
187 double x = xyLeft + i*h;
188 int idx = j*n+i;
189 mesh[idx].x = x;
190 mesh[idx].y = y;
191 if (i == 0 || i == n-1 || j == 0 || j == n-1)
192 mesh[idx].boundary = true;
193 if (i != 0)
194 mesh[idx].left = idx-1;
195 if (i != n-1)
196 mesh[idx].right = idx+1;
197 if (j != 0)
198 mesh[idx].down = idx-n;
199 if (j != n-1)
200 mesh[idx].up = idx+n;
201 if (!(eliminate_bcs && mesh[idx].boundary))
202 global_dof_indices.push_back(idx);
203 }
204 }
205
206 // Solution vector map
207 int n_global_dof = global_dof_indices.size();
208 int n_proc = comm->NumProc();
209 int proc_id = comm->MyPID();
210 int n_my_dof = n_global_dof / n_proc;
211 if (proc_id == n_proc-1)
212 n_my_dof += n_global_dof % n_proc;
213 int *my_dof = global_dof_indices.getRawPtr() + proc_id*(n_global_dof / n_proc);
214 x_map =
215 Teuchos::rcp(new Epetra_Map(n_global_dof, n_my_dof, my_dof, 0, *comm));
216
217 // Initial guess, initialized to 0.0
218 x_init = Teuchos::rcp(new Epetra_Vector(*x_map));
219 x_init->PutScalar(0.0);
220
221 // Parameter vector map
222 p_map = Teuchos::rcp(new Epetra_LocalMap(d, 0, *comm));
223
224 // Response vector map
225 g_map = Teuchos::rcp(new Epetra_LocalMap(1, 0, *comm));
226
227 // Initial parameters
228 p_init = Teuchos::rcp(new Epetra_Vector(*p_map));
229 p_init->PutScalar(0.0);
230
231 // Parameter names
232 p_names = Teuchos::rcp(new Teuchos::Array<std::string>(d));
233 for (int i=0;i<d;i++) {
234 std::stringstream ss;
235 ss << "KL Random Variable " << i+1;
236 (*p_names)[i] = ss.str();
237 }
238
239 // Build Jacobian graph
240 int NumMyElements = x_map->NumMyElements();
241 int *MyGlobalElements = x_map->MyGlobalElements();
242 graph = Teuchos::rcp(new Epetra_CrsGraph(Copy, *x_map, 5));
243 for (int i=0; i<NumMyElements; ++i ) {
244
245 // Center
246 int global_idx = MyGlobalElements[i];
247 graph->InsertGlobalIndices(global_idx, 1, &global_idx);
248
249 if (!mesh[global_idx].boundary) {
250 // Down
251 if (!(eliminate_bcs && mesh[mesh[global_idx].down].boundary))
252 graph->InsertGlobalIndices(global_idx, 1, &mesh[global_idx].down);
253
254 // Left
255 if (!(eliminate_bcs && mesh[mesh[global_idx].left].boundary))
256 graph->InsertGlobalIndices(global_idx, 1, &mesh[global_idx].left);
257
258 // Right
259 if (!(eliminate_bcs && mesh[mesh[global_idx].right].boundary))
260 graph->InsertGlobalIndices(global_idx, 1, &mesh[global_idx].right);
261
262 // Up
263 if (!(eliminate_bcs && mesh[mesh[global_idx].up].boundary))
264 graph->InsertGlobalIndices(global_idx, 1, &mesh[global_idx].up);
265 }
266 }
267 graph->FillComplete();
268 graph->OptimizeStorage();
269
270 KL_Diffusion_Func klFunc(xyLeft, xyRight, mu, s, 1.0, d);
271 if (!log_normal) {
272 // Fill coefficients of KL expansion of operator
273 if (basis == Teuchos::null) {
274 fillMatrices(klFunc, d+1);
275 }
276 else {
277 Normalized_KL_Diffusion_Func<KL_Diffusion_Func> nklFunc(klFunc, *basis);
278 fillMatrices(nklFunc, d+1);
279 }
280 }
281 else {
282 // Fill coefficients of PC expansion of operator
283 int sz = basis->size();
284 Teuchos::RCP<const Stokhos::ProductBasis<int, double> > prodbasis =
285 Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<int, double> >(
286 basis, true);
288 fillMatrices(lnFunc, sz);
289 }
290
291 // Construct deterministic operator
292 A = Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
293
294 // Construct the RHS vector.
295 b = Teuchos::rcp(new Epetra_Vector(*x_map));
296 for( int i=0 ; i<NumMyElements; ++i ) {
297 int global_idx = MyGlobalElements[i];
298 if (mesh[global_idx].boundary)
299 (*b)[i] = 0;
300 else
301 (*b)[i] = 1;
302 }
303
304 if (basis != Teuchos::null) {
305 point.resize(d);
306 basis_vals.resize(basis->size());
307 }
308}
309
310// Overridden from EpetraExt::ModelEvaluator
311
312Teuchos::RCP<const Epetra_Map>
314get_x_map() const
315{
316 return x_map;
317}
318
319Teuchos::RCP<const Epetra_Map>
321get_f_map() const
322{
323 return x_map;
324}
325
326Teuchos::RCP<const Epetra_Map>
328get_p_map(int l) const
329{
330 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
331 std::logic_error,
332 std::endl <<
333 "Error! twoD_diffusion_problem::get_p_map(): " <<
334 "Invalid parameter index l = " << l << std::endl);
335
336 return p_map;
337}
338
339Teuchos::RCP<const Epetra_Map>
341get_g_map(int l) const
342{
343 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
344 std::logic_error,
345 std::endl <<
346 "Error! twoD_diffusion_problem::get_g_map(): " <<
347 "Invalid parameter index l = " << l << std::endl);
348
349 return g_map;
350}
351
352Teuchos::RCP<const Teuchos::Array<std::string> >
354get_p_names(int l) const
355{
356 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
357 std::logic_error,
358 std::endl <<
359 "Error! twoD_diffusion_problem::get_p_names(): " <<
360 "Invalid parameter index l = " << l << std::endl);
361
362 return p_names;
363}
364
365Teuchos::RCP<const Epetra_Vector>
367get_x_init() const
368{
369 return x_init;
370}
371
372Teuchos::RCP<const Epetra_Vector>
374get_p_init(int l) const
375{
376 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
377 std::logic_error,
378 std::endl <<
379 "Error! twoD_diffusion_problem::get_p_init(): " <<
380 "Invalid parameter index l = " << l << std::endl);
381
382 return p_init;
383}
384
385Teuchos::RCP<Epetra_Operator>
387create_W() const
388{
389 Teuchos::RCP<Epetra_CrsMatrix> AA =
390 Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
391 AA->FillComplete();
392 AA->OptimizeStorage();
393 return AA;
394}
395
396void
399 const Epetra_Vector& p,
401{
402 // f = A*x - b
403 compute_A(p);
404 A->Apply(x,f);
405 f.Update(-1.0, *b, 1.0);
406}
407
408void
411 const Epetra_Vector& p,
413{
414 // J = A
415 compute_A(p);
416 Epetra_CrsMatrix& jac = dynamic_cast<Epetra_CrsMatrix&>(J);
417 jac = *A;
418 jac.FillComplete();
419 jac.OptimizeStorage();
420}
421
422void
425 const Epetra_Vector& p,
427{
428 // g = average of x
429 x.MeanValue(&g[0]);
430 g[0] *= double(x.GlobalLength()) / double(mesh.size());
431}
432
433void
439{
440 // Get stochastic expansion data
441 typedef Stokhos::Sparse3Tensor<int,double> Cijk_type;
442 Teuchos::RCP<const Cijk_type> Cijk = expn.getTripleProduct();
443 const Teuchos::Array<double>& norms = basis->norm_squared();
444
445 if (sg_kx_vec_all.size() != basis->size()) {
446 sg_kx_vec_all.resize(basis->size());
447 for (int i=0;i<basis->size();i++) {
448 sg_kx_vec_all[i] = Teuchos::rcp(new Epetra_Vector(*x_map));
449 }
450 }
451 f_sg.init(0.0);
452
453 Cijk_type::k_iterator k_begin = Cijk->k_begin();
454 Cijk_type::k_iterator k_end = Cijk->k_end();
455 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
456 int k = Stokhos::index(k_it);
457 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
458 j_it != Cijk->j_end(k_it); ++j_it) {
459 int j = Stokhos::index(j_it);
460 A_k[k]->Apply(x_sg[j],*(sg_kx_vec_all[j]));
461 }
462 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
463 j_it != Cijk->j_end(k_it); ++j_it) {
464 int j = Stokhos::index(j_it);
465 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
466 i_it != Cijk->i_end(j_it); ++i_it) {
467 int i = Stokhos::index(i_it);
468 double c = Stokhos::value(i_it); // C(i,j,k)
469 f_sg[i].Update(1.0*c/norms[i],*(sg_kx_vec_all[j]),1.0);
470 }
471 }
472 }
473 f_sg[0].Update(-1.0,*b,1.0);
474}
475
476void
481{
482 for (int i=0; i<J_sg.size(); i++) {
483 Teuchos::RCP<Epetra_CrsMatrix> jac =
484 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(J_sg.getCoeffPtr(i), true);
485 *jac = *A_k[i];
486 jac->FillComplete();
487 jac->OptimizeStorage();
488 }
489}
490
491void
496{
497 int sz = x_sg.size();
498 for (int i=0; i<sz; i++) {
499 x_sg[i].MeanValue(&g_sg[i][0]);
500 g_sg[i][0] *= double(x_sg[i].GlobalLength()) / double(mesh.size());
501 }
502}
503
504template <typename FuncT>
505void
507fillMatrices(const FuncT& func, int sz)
508{
509 int NumMyElements = x_map->NumMyElements();
510 int *MyGlobalElements = x_map->MyGlobalElements();
511 double h2 = h*h;
512 double val;
513
514 A_k.resize(sz);
515 for (int k=0; k<sz; k++) {
516 A_k[k] = Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
517 for( int i=0 ; i<NumMyElements; ++i ) {
518
519 // Center
520 int global_idx = MyGlobalElements[i];
521 if (mesh[global_idx].boundary) {
522 if (k == 0)
523 val = 1.0; //Enforce BC in mean matrix
524 else
525 val = 0.0;
526 A_k[k]->ReplaceGlobalValues(global_idx, 1, &val, &global_idx);
527 }
528 else {
529 double a_down =
530 -func(mesh[global_idx].x, mesh[global_idx].y-h/2.0, k)/h2;
531 double a_left =
532 -func(mesh[global_idx].x-h/2.0, mesh[global_idx].y, k)/h2;
533 double a_right =
534 -func(mesh[global_idx].x+h/2.0, mesh[global_idx].y, k)/h2;
535 double a_up =
536 -func(mesh[global_idx].x, mesh[global_idx].y+h/2.0, k)/h2;
537
538 // Center
539 val = -(a_down + a_left + a_right + a_up);
540 A_k[k]->ReplaceGlobalValues(global_idx, 1, &val, &global_idx);
541
542 // Down
543 if (!(eliminate_bcs && mesh[mesh[global_idx].down].boundary))
544 A_k[k]->ReplaceGlobalValues(global_idx, 1, &a_down,
545 &mesh[global_idx].down);
546
547 // Left
548 if (!(eliminate_bcs && mesh[mesh[global_idx].left].boundary))
549 A_k[k]->ReplaceGlobalValues(global_idx, 1, &a_left,
550 &mesh[global_idx].left);
551
552 // Right
553 if (!(eliminate_bcs && mesh[mesh[global_idx].right].boundary))
554 A_k[k]->ReplaceGlobalValues(global_idx, 1, &a_right,
555 &mesh[global_idx].right);
556
557 // Up
558 if (!(eliminate_bcs && mesh[mesh[global_idx].up].boundary))
559 A_k[k]->ReplaceGlobalValues(global_idx, 1, &a_up,
560 &mesh[global_idx].up);
561 }
562 }
563 A_k[k]->FillComplete();
564 A_k[k]->OptimizeStorage();
565 }
566}
567
568void
570compute_A(const Epetra_Vector& p)
571{
572 if (basis != Teuchos::null) {
573 for (int i=0; i<point.size(); i++)
574 point[i] = p[i];
575 basis->evaluateBases(point, basis_vals);
576 A->PutScalar(0.0);
577 for (int k=0;k<A_k.size();k++)
578 EpetraExt::MatrixMatrix::Add((*A_k[k]), false, basis_vals[k], *A, 1.0);
579 }
580 else {
581 *A = *(A_k[0]);
582 for (int k=1;k<A_k.size();k++)
583 EpetraExt::MatrixMatrix::Add((*A_k[k]), false, p[k-1], *A, 1.0);
584 }
585 A->FillComplete();
586 A->OptimizeStorage();
587}
Copy
expr val()
int FillComplete(bool OptimizeDataStorage=true)
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
Class representing a KL expansion of an exponential random field.
Abstract base class for multivariate orthogonal polynomials.
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const =0
Evaluate basis polynomials at given point point.
Abstract base class for orthogonal polynomial-based expansions.
virtual Teuchos::RCP< const Sparse3Tensor< ordinal_type, value_type > > getTripleProduct() const =0
Get triple product.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
Teuchos::RCP< coeff_type > getCoeffPtr(ordinal_type i)
Return ref-count pointer to coefficient i.
void init(const value_type &val)
Initialize coefficients.
ordinal_type size() const
Return size.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Teuchos::RCP< Epetra_Vector > p_init
Initial parameters.
void fillMatrices(const FuncT &func, int sz)
Fill coefficient matrix given function to evaluate diffusion coefficient.
void computeSGJacobian(const Stokhos::EpetraVectorOrthogPoly &x_sg, const Stokhos::EpetraVectorOrthogPoly &p_sg, Stokhos::EpetraOperatorOrthogPoly &J_sg)
Compute Jacobian.
Teuchos::Array< double > basis_vals
Array to store values of basis at a point.
Teuchos::RCP< Epetra_Map > x_map
Solution vector map.
Teuchos::RCP< Epetra_Vector > x_init
Initial guess.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
void compute_A(const Epetra_Vector &p)
Compute A matrix.
Teuchos::RCP< Epetra_CrsMatrix > A
Matrix to store deterministic operator.
Teuchos::RCP< Teuchos::Array< std::string > > p_names
Parameter names.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
void computeSGResponse(const Stokhos::EpetraVectorOrthogPoly &x_sg, const Stokhos::EpetraVectorOrthogPoly &p_sg, Stokhos::EpetraVectorOrthogPoly &g_sg)
Compute SG response.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
Return response function map.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::RCP< Epetra_Map > p_map
Parameter vector map.
Teuchos::Array< Teuchos::RCP< Epetra_Vector > > sg_kx_vec_all
Vectors to store matrix-vector products in SG residual calculation.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
void computeResidual(const Epetra_Vector &x, const Epetra_Vector &p, Epetra_Vector &f)
Compute residual.
Teuchos::RCP< LogNormal_Diffusion_Func< KL_Diffusion_Func > > lnFunc
twoD_diffusion_problem(const Teuchos::RCP< const Epetra_Comm > &comm, int n, int d, double s=0.1, double mu=0.2, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &basis=Teuchos::null, bool log_normal=false, bool eliminate_bcs=false)
Constructor.
Teuchos::Array< MeshPoint > mesh
Teuchos::Array< double > point
Array to store a point for basis evaluation.
void computeJacobian(const Epetra_Vector &x, const Epetra_Vector &p, Epetra_Operator &J)
Compute Jacobian.
Teuchos::Array< Teuchos::RCP< Epetra_CrsMatrix > > A_k
KL coefficients of operator
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
Teuchos::RCP< KL_Diffusion_Func > klFunc
Teuchos::RCP< Epetra_Vector > b
Deterministic RHS.
Teuchos::RCP< Epetra_CrsGraph > graph
Jacobian graph.
void computeResponse(const Epetra_Vector &x, const Epetra_Vector &p, Epetra_Vector &g)
Compute response.
Teuchos::RCP< Epetra_Map > g_map
Response vector map.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
void computeSGResidual(const Stokhos::EpetraVectorOrthogPoly &x_sg, const Stokhos::EpetraVectorOrthogPoly &p_sg, Stokhos::OrthogPolyExpansion< int, double > &expn, Stokhos::EpetraVectorOrthogPoly &f_sg)
Compute SG residual.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:265
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition: csr_vector.h:267