Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
twoD_diffusion_problem_tpetra.hpp
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 #ifndef TWOD_DIFFUSION_PROBLEM_TPETRA_HPP
43 #define TWOD_DIFFUSION_PROBLEM_TPETRA_HPP
44 
45 #include "Teuchos_RCP.hpp"
46 #include "Teuchos_Array.hpp"
47 #include "Teuchos_Comm.hpp"
48 
49 #include "Tpetra_DefaultPlatform.hpp"
50 #include "Tpetra_Map.hpp"
51 #include "Tpetra_Vector.hpp"
52 #include "Tpetra_Import.hpp"
53 #include "Tpetra_CrsGraph.hpp"
54 #include "Tpetra_CrsMatrix.hpp"
55 
56 #include "Stokhos.hpp"
58 //#include "Stokhos_Epetra.hpp"
59 
61 template <typename Scalar,
62  typename MeshScalar,
63  typename BasisScalar,
64  typename LocalOrdinal,
65  typename GlobalOrdinal,
66  typename Node>
68 public:
69 
70  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
72  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Operator;
73  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
74  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
75  typedef Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Import;
76 
79  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
81  BasisScalar s = 0.1, BasisScalar mu = 0.2,
82  bool log_normal = false,
83  bool eliminate_bcs = false);
84 
87 
89  Teuchos::RCP<const Tpetra_Map> get_x_map() const;
90 
92  Teuchos::RCP<const Tpetra_Map> get_f_map() const;
93 
95  Teuchos::RCP<const Tpetra_Map> get_p_map(LocalOrdinal l) const;
96 
98  Teuchos::RCP<const Tpetra_Map> get_g_map(LocalOrdinal j) const;
99 
101  Teuchos::RCP<const Teuchos::Array<std::string> >
102  get_p_names(LocalOrdinal l) const;
103 
105  Teuchos::RCP<const Tpetra_Vector> get_x_init() const;
106 
108  Teuchos::RCP<const Tpetra_Vector> get_p_init(LocalOrdinal l) const;
109 
111  Teuchos::RCP<Tpetra_CrsMatrix> create_W() const;
112 
114  void computeResidual(const Tpetra_Vector& x,
115  const Tpetra_Vector& p,
116  Tpetra_Vector& f);
117 
119  void computeJacobian(const Tpetra_Vector& x,
120  const Tpetra_Vector& p,
121  Tpetra_CrsMatrix& J);
122 
124  void computeResponse(const Tpetra_Vector& x,
125  const Tpetra_Vector& p,
126  Tpetra_Vector& g);
127 
129 
130 protected:
131 
133  template <typename FuncT>
134  void computeA(const FuncT& func, const Tpetra_Vector& p,
135  Tpetra_CrsMatrix& jac);
136 
137 protected:
138 
139  // Mesh
140  MeshScalar h;
141  struct MeshPoint {
142  MeshPoint() : up(-1), down(-1), left(-1), right(-1), boundary(false) {}
143  MeshScalar x, y;
145  bool boundary;
146  };
147  Teuchos::Array<MeshPoint> mesh;
148  Teuchos::Array<GlobalOrdinal> bcIndices;
149  bool log_normal;
150  bool eliminate_bcs;
151 
153  Teuchos::RCP<const Tpetra_Map> x_map;
154 
156  Teuchos::RCP<Tpetra_Import> importer;
157 
159  Teuchos::RCP<Tpetra_Vector> x_init;
160 
162  Teuchos::RCP<const Tpetra_Map> p_map;
163 
165  Teuchos::RCP<const Tpetra_Map> g_map;
166 
168  Teuchos::RCP<Tpetra_Vector> p_init;
169 
171  Teuchos::RCP< Teuchos::Array<std::string> > p_names;
172 
174  Teuchos::RCP<Tpetra_CrsGraph> graph;
175 
177  Teuchos::RCP<Tpetra_Vector> b;
178 
180  Teuchos::RCP<Tpetra_CrsMatrix> A;
181 
182  // Functor representing a diffusion function given by a KL expansion
184  mutable Teuchos::Array<MeshScalar> point;
185  Teuchos::RCP< Stokhos::KL::ExponentialRandomField<MeshScalar> > rf;
186  KL_Diffusion_Func(MeshScalar xyLeft, MeshScalar xyRight,
187  BasisScalar mean, BasisScalar std_dev,
188  MeshScalar L, LocalOrdinal num_KL);
189  Scalar operator() (MeshScalar x, MeshScalar y,
190  const Teuchos::Array<Scalar>& rv) const;
191  };
192 
193  // Functor representing a diffusion function given by a log-normal
194  // PC expansion
195  template <typename DiffusionFunc>
197  const DiffusionFunc& func;
198  LogNormal_Diffusion_Func(const DiffusionFunc& func_) : func(func_) {}
199  Scalar operator() (MeshScalar x, MeshScalar y,
200  const Teuchos::Array<Scalar>& rv) const {
201  return std::exp(func(x,y,rv));
202  }
203  };
204 
205  Teuchos::RCP<KL_Diffusion_Func> klFunc;
206  Teuchos::RCP<LogNormal_Diffusion_Func<KL_Diffusion_Func> > lnFunc;
207 
208 };
209 
211 
212 #endif // TWOD_DIFFUSION_ME_HPP
Scalar operator()(MeshScalar x, MeshScalar y, const Teuchos::Array< Scalar > &rv) const
Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > Tpetra_Map
KL_Diffusion_Func(MeshScalar xyLeft, MeshScalar xyRight, BasisScalar mean, BasisScalar std_dev, MeshScalar L, LocalOrdinal num_KL)
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > Tpetra_Operator
Teuchos::RCP< const Tpetra_Map > p_map
Parameter vector map.
void computeResidual(const Epetra_Vector &x, const Epetra_Vector &p, Epetra_Vector &f)
Compute residual.
void computeResponse(const Epetra_Vector &x, const Epetra_Vector &p, Epetra_Vector &g)
Compute response.
Teuchos::RCP< LogNormal_Diffusion_Func< KL_Diffusion_Func > > lnFunc
void computeJacobian(const Epetra_Vector &x, const Epetra_Vector &p, Epetra_Operator &J)
Compute Jacobian.
Teuchos::RCP< Tpetra_Vector > x_init
Initial guess.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< const Tpetra_Map > g_map
Response vector map.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
Tpetra::Import< LocalOrdinal, GlobalOrdinal, Node > Tpetra_Import
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
Teuchos::Array< MeshPoint > mesh
Teuchos::RCP< Stokhos::KL::ExponentialRandomField< MeshScalar > > rf
Teuchos::RCP< KL_Diffusion_Func > klFunc
Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > Tpetra_CrsGraph
KokkosClassic::DefaultNode::DefaultNodeType Node
Teuchos::RCP< const Tpetra_Map > x_map
Solution vector map.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
expr expr expr expr j
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.
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
Teuchos::Array< GlobalOrdinal > bcIndices
Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Tpetra_Vector
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
Return response function map.
Teuchos::RCP< Tpetra_CrsMatrix > A
Matrix to store operator.
Scalar operator()(MeshScalar x, MeshScalar y, const Teuchos::Array< Scalar > &rv) const
Teuchos::RCP< Tpetra_Import > importer
Importer to overlapped distribution.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Teuchos::RCP< Tpetra_Vector > b
RHS.
pce_type Scalar
Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Tpetra_CrsMatrix
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
Teuchos::RCP< Teuchos::Array< std::string > > p_names
Parameter names.
Teuchos::RCP< Tpetra_Vector > p_init
Initial parameters.
void computeA(const FuncT &func, const Tpetra_Vector &p, Tpetra_CrsMatrix &jac)
Fill coefficient matrix given function to evaluate diffusion coefficient.
Teuchos::RCP< Tpetra_CrsGraph > graph
Jacobian graph.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition: csr_vector.h:267
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
A linear 2-D diffusion problem.