Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
twoD_diffusion_ME.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
42#include "twoD_diffusion_ME.hpp"
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, Kokkos::DefaultHostExecutionSpace> > 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 =
79 }
80
81 ~KL_Diffusion_Func() {
82 }
83
84 double operator() (double x, double y, int k) const {
85 if (k == 0)
86 return mean;
87 point[0] = x;
88 point[1] = y;
89 return rf->evaluate_eigenfunction(point, k-1);
90 }
91 };
92
93 // Functor representing a diffusion function given by a KL expansion
94 // where the basis is normalized
95 template <typename DiffusionFunc>
96 struct Normalized_KL_Diffusion_Func {
97 const DiffusionFunc& func;
98 int d;
99 Teuchos::Array<double> psi_0, psi_1;
100
101 Normalized_KL_Diffusion_Func(
102 const DiffusionFunc& func_,
104 func(func_),
105 d(basis.dimension()),
106 psi_0(basis.size()),
107 psi_1(basis.size())
108 {
109 Teuchos::Array<double> zero(d), one(d);
110 for(int k=0; k<d; k++) {
111 zero[k] = 0.0;
112 one[k] = 1.0;
113 }
114 basis.evaluateBases(zero, psi_0);
115 basis.evaluateBases(one, psi_1);
116 }
117
118 double operator() (double x, double y, int k) const {
119 if (k == 0) {
120 double val = func(x, y, 0);
121 for (int i=1; i<=d; i++)
122 val -= psi_0[i]/(psi_1[i]-psi_0[i])*func(x, y, i);
123 val /= psi_0[0];
124 return val;
125 }
126 else
127 return 1.0/(psi_1[k]-psi_0[k])*func(x, y, k);
128 }
129 };
130
131 // Functor representing a diffusion function given by a log-normal PC expansion
132 template <typename DiffusionFunc>
133 struct LogNormal_Diffusion_Func {
134 double mean;
135 const DiffusionFunc& func;
136 Teuchos::RCP<const Stokhos::ProductBasis<int, double> > prodbasis;
137
138 LogNormal_Diffusion_Func(
139 double mean_,
140 const DiffusionFunc& func_,
141 const Teuchos::RCP<const Stokhos::ProductBasis<int, double> > prodbasis_)
142 : mean(mean_), func(func_), prodbasis(prodbasis_) {}
143
144 double operator() (double x, double y, int k) const {
145 int d = prodbasis->dimension();
146 const Teuchos::Array<double>& norms = prodbasis->norm_squared();
147 const Stokhos::MultiIndex<int> multiIndex = prodbasis->term(k);
148 double sum_g = 0.0, efval;
149 for (int l=0; l<d; l++) {
150 sum_g += std::pow(func(x,y,l+1),2);
151 }
152 efval = std::exp(mean + 0.5*sum_g)/norms[k];
153 for (int l=0; l<d; l++) {
154 efval *= std::pow(func(x,y,l+1),multiIndex[l]);
155 }
156 return efval;
157 }
158 };
159
160}
161
164 const Teuchos::RCP<const Epetra_Comm>& comm, int n, int d,
165 double s, double mu,
166 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& basis_,
167 bool log_normal_,
168 bool eliminate_bcs_,
169 const Teuchos::RCP<Teuchos::ParameterList>& precParams_) :
170 mesh(n*n),
171 basis(basis_),
172 log_normal(log_normal_),
173 eliminate_bcs(eliminate_bcs_),
174 precParams(precParams_)
175{
176 Kokkos::initialize();
177
179 // Construct the mesh.
180 // The mesh is uniform and the nodes are numbered
181 // LEFT to RIGHT, DOWN to UP.
182 //
183 // 5-6-7-8-9
184 // | | | | |
185 // 0-1-2-3-4
187 double xyLeft = -.5;
188 double xyRight = .5;
189 h = (xyRight - xyLeft)/((double)(n-1));
190 Teuchos::Array<int> global_dof_indices;
191 for (int j=0; j<n; j++) {
192 double y = xyLeft + j*h;
193 for (int i=0; i<n; i++) {
194 double x = xyLeft + i*h;
195 int idx = j*n+i;
196 mesh[idx].x = x;
197 mesh[idx].y = y;
198 if (i == 0 || i == n-1 || j == 0 || j == n-1)
199 mesh[idx].boundary = true;
200 if (i != 0)
201 mesh[idx].left = idx-1;
202 if (i != n-1)
203 mesh[idx].right = idx+1;
204 if (j != 0)
205 mesh[idx].down = idx-n;
206 if (j != n-1)
207 mesh[idx].up = idx+n;
208 if (!(eliminate_bcs && mesh[idx].boundary))
209 global_dof_indices.push_back(idx);
210 }
211 }
212
213 // Solution vector map
214 int n_global_dof = global_dof_indices.size();
215 int n_proc = comm->NumProc();
216 int proc_id = comm->MyPID();
217 int n_my_dof = n_global_dof / n_proc;
218 if (proc_id == n_proc-1)
219 n_my_dof += n_global_dof % n_proc;
220 int *my_dof = global_dof_indices.getRawPtr() + proc_id*(n_global_dof / n_proc);
221 x_map =
222 Teuchos::rcp(new Epetra_Map(n_global_dof, n_my_dof, my_dof, 0, *comm));
223
224 // Initial guess, initialized to 0.0
225 x_init = Teuchos::rcp(new Epetra_Vector(*x_map));
226 x_init->PutScalar(0.0);
227
228 // Parameter vector map
229 p_map = Teuchos::rcp(new Epetra_LocalMap(d, 0, *comm));
230
231 // Response vector map
232 g_map = Teuchos::rcp(new Epetra_LocalMap(1, 0, *comm));
233
234 // Initial parameters
235 p_init = Teuchos::rcp(new Epetra_Vector(*p_map));
236 p_init->PutScalar(0.0);
237
238 // Parameter names
239 p_names = Teuchos::rcp(new Teuchos::Array<std::string>(d));
240 for (int i=0;i<d;i++) {
241 std::stringstream ss;
242 ss << "KL Random Variable " << i+1;
243 (*p_names)[i] = ss.str();
244 }
245
246 // Build Jacobian graph
247 int NumMyElements = x_map->NumMyElements();
248 int *MyGlobalElements = x_map->MyGlobalElements();
249 graph = Teuchos::rcp(new Epetra_CrsGraph(Copy, *x_map, 5));
250 for (int i=0; i<NumMyElements; ++i ) {
251
252 // Center
253 int global_idx = MyGlobalElements[i];
254 graph->InsertGlobalIndices(global_idx, 1, &global_idx);
255
256 if (!mesh[global_idx].boundary) {
257 // Down
258 if (!(eliminate_bcs && mesh[mesh[global_idx].down].boundary))
259 graph->InsertGlobalIndices(global_idx, 1, &mesh[global_idx].down);
260
261 // Left
262 if (!(eliminate_bcs && mesh[mesh[global_idx].left].boundary))
263 graph->InsertGlobalIndices(global_idx, 1, &mesh[global_idx].left);
264
265 // Right
266 if (!(eliminate_bcs && mesh[mesh[global_idx].right].boundary))
267 graph->InsertGlobalIndices(global_idx, 1, &mesh[global_idx].right);
268
269 // Up
270 if (!(eliminate_bcs && mesh[mesh[global_idx].up].boundary))
271 graph->InsertGlobalIndices(global_idx, 1, &mesh[global_idx].up);
272 }
273 }
274 graph->FillComplete();
275 graph->OptimizeStorage();
276
277 KL_Diffusion_Func klFunc(xyLeft, xyRight, mu, s, 1.0, d);
278 if (!log_normal) {
279 // Fill coefficients of KL expansion of operator
280 if (basis == Teuchos::null) {
281 fillMatrices(klFunc, d+1);
282 }
283 else {
284 Normalized_KL_Diffusion_Func<KL_Diffusion_Func> nklFunc(klFunc, *basis);
285 fillMatrices(nklFunc, d+1);
286 }
287 }
288 else {
289 // Fill coefficients of PC expansion of operator
290 int sz = basis->size();
291 Teuchos::RCP<const Stokhos::ProductBasis<int, double> > prodbasis =
292 Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<int, double> >(
293 basis, true);
294 LogNormal_Diffusion_Func<KL_Diffusion_Func> lnFunc(mu, klFunc, prodbasis);
295 fillMatrices(lnFunc, sz);
296 }
297
298 // Construct deterministic operator
299 A = Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
300
301 // Construct the RHS vector.
302 b = Teuchos::rcp(new Epetra_Vector(*x_map));
303 for( int i=0 ; i<NumMyElements; ++i ) {
304 int global_idx = MyGlobalElements[i];
305 if (mesh[global_idx].boundary)
306 (*b)[i] = 0;
307 else
308 (*b)[i] = 1;
309 }
310
311 if (basis != Teuchos::null) {
312 point.resize(d);
313 basis_vals.resize(basis->size());
314 }
315
316 if (precParams != Teuchos::null) {
317 std::string name = precParams->get("Preconditioner Type", "Ifpack");
318 Teuchos::RCP<Teuchos::ParameterList> p =
319 Teuchos::rcp(&(precParams->sublist("Preconditioner Parameters")), false);
321 Teuchos::rcp(new Stokhos::PreconditionerFactory(name, p));
322 }
323}
324
327{
328 Kokkos::finalize();
329}
330
331// Overridden from EpetraExt::ModelEvaluator
332
333Teuchos::RCP<const Epetra_Map>
335get_x_map() const
336{
337 return x_map;
338}
339
340Teuchos::RCP<const Epetra_Map>
342get_f_map() const
343{
344 return x_map;
345}
346
347Teuchos::RCP<const Epetra_Map>
349get_p_map(int l) const
350{
351 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
352 std::logic_error,
353 std::endl <<
354 "Error! twoD_diffusion_ME::get_p_map(): " <<
355 "Invalid parameter index l = " << l << std::endl);
356
357 return p_map;
358}
359
360Teuchos::RCP<const Epetra_Map>
362get_g_map(int l) const
363{
364 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
365 std::logic_error,
366 std::endl <<
367 "Error! twoD_diffusion_ME::get_g_map(): " <<
368 "Invalid parameter index l = " << l << std::endl);
369
370 return g_map;
371}
372
373Teuchos::RCP<const Teuchos::Array<std::string> >
375get_p_names(int l) const
376{
377 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
378 std::logic_error,
379 std::endl <<
380 "Error! twoD_diffusion_ME::get_p_names(): " <<
381 "Invalid parameter index l = " << l << std::endl);
382
383 return p_names;
384}
385
386Teuchos::RCP<const Epetra_Vector>
388get_x_init() const
389{
390 return x_init;
391}
392
393Teuchos::RCP<const Epetra_Vector>
395get_p_init(int l) const
396{
397 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
398 std::logic_error,
399 std::endl <<
400 "Error! twoD_diffusion_ME::get_p_init(): " <<
401 "Invalid parameter index l = " << l << std::endl);
402
403 return p_init;
404}
405
406Teuchos::RCP<Epetra_Operator>
408create_W() const
409{
410 Teuchos::RCP<Epetra_CrsMatrix> AA =
411 Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
412 AA->FillComplete();
413 AA->OptimizeStorage();
414 return AA;
415}
416
417Teuchos::RCP<EpetraExt::ModelEvaluator::Preconditioner>
419create_WPrec() const
420{
421 if (precFactory != Teuchos::null) {
422 Teuchos::RCP<Epetra_Operator> precOp =
423 precFactory->compute(A, false);
424 return Teuchos::rcp(new EpetraExt::ModelEvaluator::Preconditioner(precOp,
425 true));
426 }
427 return Teuchos::null;
428}
429
430EpetraExt::ModelEvaluator::InArgs
432createInArgs() const
433{
434 InArgsSetup inArgs;
435 inArgs.setModelEvalDescription("TwoD Diffusion Model Evaluator");
436
437 // Deterministic InArgs
438 inArgs.setSupports(IN_ARG_x,true);
439 inArgs.set_Np(1); // 1 parameter vector
440
441 // Stochastic InArgs
442 inArgs.setSupports(IN_ARG_x_sg,true);
443 inArgs.setSupports(IN_ARG_p_sg, 0, true); // 1 SG parameter vector
444 inArgs.setSupports(IN_ARG_sg_basis,true);
445 inArgs.setSupports(IN_ARG_sg_quadrature,true);
446 inArgs.setSupports(IN_ARG_sg_expansion,true);
447
448 // Multipoint InArgs
449 inArgs.setSupports(IN_ARG_x_mp,true);
450 inArgs.setSupports(IN_ARG_p_mp, 0, true); // 1 MP parameter vector
451
452 return inArgs;
453}
454
455EpetraExt::ModelEvaluator::OutArgs
457createOutArgs() const
458{
459 OutArgsSetup outArgs;
460 outArgs.setModelEvalDescription("TwoD Diffusion Model Evaluator");
461
462 // Deterministic OutArgs
463 outArgs.set_Np_Ng(1, 1);
464 outArgs.setSupports(OUT_ARG_f,true);
465 outArgs.setSupports(OUT_ARG_W,true);
466 if (precFactory != Teuchos::null)
467 outArgs.setSupports(OUT_ARG_WPrec,true);
468
469 // Stochastic OutArgs
470 outArgs.setSupports(OUT_ARG_f_sg,true);
471 outArgs.setSupports(OUT_ARG_W_sg,true);
472 outArgs.setSupports(OUT_ARG_g_sg, 0, true);
473
474 // Multipoint OutArgs
475 outArgs.setSupports(OUT_ARG_f_mp,true);
476 outArgs.setSupports(OUT_ARG_W_mp,true);
477 outArgs.setSupports(OUT_ARG_g_mp, 0, true);
478
479 return outArgs;
480}
481
482void
484evalModel(const InArgs& inArgs, const OutArgs& outArgs) const
485{
486
487 //
488 // Determinisic calculation
489 //
490
491 // Solution vector
492 Teuchos::RCP<const Epetra_Vector> det_x = inArgs.get_x();
493
494 // Parameters
495 Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(0);
496 if (p == Teuchos::null)
497 p = p_init;
498
499 Teuchos::RCP<Epetra_Vector> f = outArgs.get_f();
500 Teuchos::RCP<Epetra_Operator> W = outArgs.get_W();
501 Teuchos::RCP<Epetra_Operator> WPrec = outArgs.get_WPrec();
502 if (f != Teuchos::null || W != Teuchos::null || WPrec != Teuchos::null) {
503 if (basis != Teuchos::null) {
504 for (int i=0; i<point.size(); i++)
505 point[i] = (*p)[i];
506 basis->evaluateBases(point, basis_vals);
507 A->PutScalar(0.0);
508 for (int k=0;k<A_k.size();k++)
509 EpetraExt::MatrixMatrix::Add((*A_k[k]), false, basis_vals[k], *A, 1.0);
510 }
511 else {
512 *A = *(A_k[0]);
513 for (int k=1;k<A_k.size();k++)
514 EpetraExt::MatrixMatrix::Add((*A_k[k]), false, (*p)[k-1], *A, 1.0);
515 }
516 A->FillComplete();
517 A->OptimizeStorage();
518 }
519
520 // Residual
521 if (f != Teuchos::null) {
522 Teuchos::RCP<Epetra_Vector> kx = Teuchos::rcp(new Epetra_Vector(*x_map));
523 A->Apply(*det_x,*kx);
524 f->Update(1.0,*kx,-1.0, *b, 0.0);
525 }
526
527 // Jacobian
528 if (W != Teuchos::null) {
529 Teuchos::RCP<Epetra_CrsMatrix> jac =
530 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W, true);
531 *jac = *A;
532 jac->FillComplete();
533 jac->OptimizeStorage();
534 }
535
536 // Preconditioner
537 if (WPrec != Teuchos::null)
538 precFactory->recompute(A, WPrec);
539
540 // Responses (mean value)
541 Teuchos::RCP<Epetra_Vector> g = outArgs.get_g(0);
542 if (g != Teuchos::null) {
543 (det_x->MeanValue(&(*g)[0]));
544 (*g)[0] *= double(det_x->GlobalLength()) / double(mesh.size());
545 }
546
547 //
548 // Stochastic Galerkin calculation
549 //
550
551 // Stochastic solution vector
552 InArgs::sg_const_vector_t x_sg = inArgs.get_x_sg();
553
554 // Stochastic parameters
555 InArgs::sg_const_vector_t p_sg = inArgs.get_p_sg(0);
556
557 // Stochastic residual
558 OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
559 if (f_sg != Teuchos::null) {
560
561 // Get stochastic expansion data
562 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expn =
563 inArgs.get_sg_expansion();
564 typedef Stokhos::Sparse3Tensor<int,double> Cijk_type;
565 Teuchos::RCP<const Cijk_type> Cijk = expn->getTripleProduct();
566 const Teuchos::Array<double>& norms = basis->norm_squared();
567
568 if (sg_kx_vec_all.size() != basis->size()) {
569 sg_kx_vec_all.resize(basis->size());
570 for (int i=0;i<basis->size();i++) {
571 sg_kx_vec_all[i] = Teuchos::rcp(new Epetra_Vector(*x_map));
572 }
573 }
574 f_sg->init(0.0);
575
576 Cijk_type::k_iterator k_begin = Cijk->k_begin();
577 Cijk_type::k_iterator k_end = Cijk->k_end();
578 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
579 int k = Stokhos::index(k_it);
580 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
581 j_it != Cijk->j_end(k_it); ++j_it) {
582 int j = Stokhos::index(j_it);
583 A_k[k]->Apply((*x_sg)[j],*(sg_kx_vec_all[j]));
584 }
585 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
586 j_it != Cijk->j_end(k_it); ++j_it) {
587 int j = Stokhos::index(j_it);
588 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
589 i_it != Cijk->i_end(j_it); ++i_it) {
590 int i = Stokhos::index(i_it);
591 double c = Stokhos::value(i_it); // C(i,j,k)
592 (*f_sg)[i].Update(1.0*c/norms[i],*(sg_kx_vec_all[j]),1.0);
593 }
594 }
595 } //End
596 (*f_sg)[0].Update(-1.0,*b,1.0);
597 }
598
599 // Stochastic Jacobian
600 OutArgs::sg_operator_t W_sg = outArgs.get_W_sg();
601 if (W_sg != Teuchos::null) {
602 for (int i=0; i<W_sg->size(); i++) {
603 Teuchos::RCP<Epetra_CrsMatrix> jac =
604 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_sg->getCoeffPtr(i), true);
605 *jac = *A_k[i];
606 jac->FillComplete();
607 jac->OptimizeStorage();
608 }
609 }
610
611 // Stochastic responses
612 Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > g_sg =
613 outArgs.get_g_sg(0);
614 if (g_sg != Teuchos::null) {
615 int sz = x_sg->size();
616 for (int i=0; i<sz; i++) {
617 (*x_sg)[i].MeanValue(&(*g_sg)[i][0]);
618 (*g_sg)[i][0] *= double((*x_sg)[i].GlobalLength()) / double(mesh.size());
619 }
620 }
621
622 //
623 // Multi-point calculation
624 //
625
626 // Stochastic solution vector
627 mp_const_vector_t x_mp = inArgs.get_x_mp();
628
629 // Stochastic parameters
630 mp_const_vector_t p_mp = inArgs.get_p_mp(0);
631
632 // Stochastic residual
633 mp_vector_t f_mp = outArgs.get_f_mp();
634 mp_operator_t W_mp = outArgs.get_W_mp();
635 if (f_mp != Teuchos::null || W_mp != Teuchos::null) {
636 int num_mp = x_mp->size();
637 for (int i=0; i<num_mp; i++) {
638 // Compute operator
639 if (basis != Teuchos::null) {
640 for (int k=0; k<point.size(); k++)
641 point[k] = (*p_mp)[i][k];
642 basis->evaluateBases(point, basis_vals);
643 A->PutScalar(0.0);
644 for (int k=0;k<A_k.size();k++)
645 EpetraExt::MatrixMatrix::Add((*A_k[k]), false, basis_vals[k], *A,
646 1.0);
647 }
648 else {
649 *A = *(A_k[0]);
650 for (int k=1;k<A_k.size();k++)
651 EpetraExt::MatrixMatrix::Add((*A_k[k]), false, (*p_mp)[i][k-1], *A,
652 1.0);
653 }
654
655 A->FillComplete();
656 A->OptimizeStorage();
657
658 // Compute residual
659 if (f_mp != Teuchos::null) {
660 A->Apply((*x_mp)[i], (*f_mp)[i]);
661 (*f_mp)[i].Update(-1.0, *b, 1.0);
662 }
663
664 // Copy operator
665 if (W_mp != Teuchos::null) {
666 Teuchos::RCP<Epetra_CrsMatrix> jac =
667 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_mp->getCoeffPtr(i),
668 true);
669 *jac = *A;
670 jac->FillComplete();
671 jac->OptimizeStorage();
672 }
673 }
674 }
675
676 // Multipoint responses
677 mp_vector_t g_mp = outArgs.get_g_mp(0);
678 if (g_mp != Teuchos::null) {
679 int sz = x_mp->size();
680 for (int i=0; i<sz; i++) {
681 (*x_mp)[i].MeanValue(&(*g_mp)[i][0]);
682 (*g_mp)[i][0] *= double((*x_mp)[i].GlobalLength()) / double(mesh.size());
683 }
684 }
685
686}
687
Copy
expr val()
Class representing a KL expansion of an exponential random field.
A multidimensional index.
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.
An class for building preconditioners.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Teuchos::RCP< Epetra_Vector > p_init
Initial parameters.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
InArgs createInArgs() const
Create InArgs.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
Teuchos::Array< MeshPoint > mesh
Teuchos::Array< double > point
Array to store a point for basis evaluation.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
OutArgs createOutArgs() const
Create OutArgs.
twoD_diffusion_ME(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, const Teuchos::RCP< Teuchos::ParameterList > &precParams=Teuchos::null)
Constructor.
Teuchos::RCP< Epetra_CrsMatrix > A
Matrix to store deterministic operator.
Teuchos::Array< Teuchos::RCP< Epetra_CrsMatrix > > A_k
KL coefficients of operator
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::RCP< Epetra_Map > x_map
Solution vector map.
~twoD_diffusion_ME()
Destructor.
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
Return response function map.
Teuchos::RCP< Epetra_Map > g_map
Response vector map.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::RCP< Epetra_Map > p_map
Parameter vector map.
Teuchos::RCP< Epetra_Vector > x_init
Initial guess.
Teuchos::Array< Teuchos::RCP< Epetra_Vector > > sg_kx_vec_all
Vectors to store matrix-vector products in SG residual calculation.
Teuchos::RCP< Stokhos::AbstractPreconditionerFactory > precFactory
Teuchos::RCP< Teuchos::ParameterList > precParams
Teuchos::RCP< Epetra_CrsGraph > graph
Jacobian graph.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::RCP< EpetraExt::ModelEvaluator::Preconditioner > create_WPrec() const
Create preconditioner for W.
void fillMatrices(const FuncT &func, int sz)
Fill coefficient matrix given function to evaluate diffusion coefficient.
Teuchos::RCP< Teuchos::Array< std::string > > p_names
Parameter names.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::Array< double > basis_vals
Array to store values of basis at a point.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
Teuchos::RCP< Epetra_Vector > b
Deterministic RHS.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
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