Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ExplicitModelEvaluator_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Panzer: A partial differential equation assembly
5// engine for strongly coupled complex multiphysics systems
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39// Eric C. Cyr (eccyr@sandia.gov)
40// ***********************************************************************
41// @HEADER
42
43#ifndef PANZER_EXPLICIT_MODEL_EVALUATOR_IMPL_HPP
44#define PANZER_EXPLICIT_MODEL_EVALUATOR_IMPL_HPP
45
46#include "Thyra_DefaultDiagonalLinearOp.hpp"
47#include "Thyra_LinearOpWithSolveFactoryBase.hpp"
48
49#include "PanzerDiscFE_config.hpp"
50#ifdef PANZER_HAVE_EPETRA_STACK
51#include "Thyra_EpetraModelEvaluator.hpp"
52#include "Thyra_get_Epetra_Operator.hpp"
54#endif
55
56namespace panzer {
57
58template<typename Scalar>
60ExplicitModelEvaluator(const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > & model,
61 bool constantMassMatrix,
62 bool useLumpedMass,
63 bool applyMassInverse)
64 : Thyra::ModelEvaluatorDelegatorBase<Scalar>(model)
65 , constantMassMatrix_(constantMassMatrix)
66 , massLumping_(useLumpedMass)
67{
68 using Teuchos::RCP;
69 using Teuchos::rcp_dynamic_cast;
70
71 this->applyMassInverse_ = applyMassInverse;
72
73 // extract a panzer::ModelEvaluator if appropriate
74 panzerModel_ = rcp_dynamic_cast<panzer::ModelEvaluator<Scalar> >(model);
75
76#ifdef PANZER_HAVE_EPETRA_STACK
77 // extract a panzer::ModelEvaluator_Epetra if appropriate
78 RCP<Thyra::EpetraModelEvaluator> epME = rcp_dynamic_cast<Thyra::EpetraModelEvaluator>(model);
79 if(epME!=Teuchos::null)
80 panzerEpetraModel_ = rcp_dynamic_cast<const panzer::ModelEvaluator_Epetra>(epME->getEpetraModel());
81#endif
82 // note at this point its possible that panzerModel_ = panzerEpetraModel_ = Teuchos::null
83
85}
86
87template<typename Scalar>
88Thyra::ModelEvaluatorBase::InArgs<Scalar> ExplicitModelEvaluator<Scalar>::
89getNominalValues() const
90{
91 typedef Thyra::ModelEvaluatorBase MEB;
92
93 MEB::InArgs<Scalar> nomVals = createInArgs();
94 nomVals.setArgs(this->getUnderlyingModel()->getNominalValues(),true);
95
96 return nomVals;
97}
98
99template<typename Scalar>
100Thyra::ModelEvaluatorBase::InArgs<Scalar> ExplicitModelEvaluator<Scalar>::
101createInArgs() const
102{
103 return prototypeInArgs_;
104}
105
106template<typename Scalar>
107Thyra::ModelEvaluatorBase::OutArgs<Scalar> ExplicitModelEvaluator<Scalar>::
108createOutArgs() const
109{
110 return prototypeOutArgs_;
111}
112
113template<typename Scalar>
114Teuchos::RCP<panzer::ModelEvaluator<Scalar> > ExplicitModelEvaluator<Scalar>::
116{
117 return Teuchos::rcp_dynamic_cast<panzer::ModelEvaluator<Scalar> >(this->getNonconstUnderlyingModel());
118}
119
120template<typename Scalar>
122evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
123 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
124{
125 typedef Thyra::ModelEvaluatorBase MEB;
126 using Teuchos::RCP;
127 RCP<const Thyra::ModelEvaluator<Scalar> > under_me = this->getUnderlyingModel();
128
129 MEB::InArgs<Scalar> under_inArgs = under_me->createInArgs();
130 under_inArgs.setArgs(inArgs);
131
132 // read in the supplied time derivative in case it is needed by the explicit model (e.g. for stabilization) and make sure alpha is set to zero
133 under_inArgs.set_x_dot(inArgs.get_x_dot());
134 under_inArgs.set_alpha(0.0);
135
136 Teuchos::RCP<Thyra::VectorBase<Scalar> > f = outArgs.get_f();
137
138 MEB::OutArgs<Scalar> under_outArgs = under_me->createOutArgs();
139 under_outArgs.setArgs(outArgs);
140 if(f!=Teuchos::null) {
141 // build a scrap vector that will contain the weak residual
142 if(scrap_f_==Teuchos::null)
143 scrap_f_ = Thyra::createMember(*under_me->get_f_space());
144
145 Thyra::assign(scrap_f_.ptr(),0.0);
146 under_outArgs.set_f(scrap_f_);
147 }
148
149 under_me->evalModel(under_inArgs,under_outArgs);
150
151 // build the mass matrix
152 if(invMassMatrix_==Teuchos::null || constantMassMatrix_==false)
153 buildInverseMassMatrix(inArgs);
154
155 if(f!=Teuchos::null)
156 Thyra::Vt_S(scrap_f_.ptr(),-1.0);
157
158 // invert the mass matrix
159 if(f!=Teuchos::null && this->applyMassInverse_) {
160 Thyra::apply(*invMassMatrix_,Thyra::NOTRANS,*scrap_f_,f.ptr());
161 }
162 else if(f!=Teuchos::null){
163 Thyra::V_V(f.ptr(),*scrap_f_);
164 }
165}
166
167template<typename Scalar>
169buildInverseMassMatrix(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs) const
170{
171 typedef Thyra::ModelEvaluatorBase MEB;
172 using Teuchos::RCP;
173 using Thyra::createMember;
174
175 RCP<const Thyra::ModelEvaluator<Scalar> > me = this->getUnderlyingModel();
176
177 // first allocate space for the mass matrix
178 mass_ = me->create_W_op();
179
180 // intialize a zero to get rid of the x-dot
181 if(zero_==Teuchos::null) {
182 zero_ = Thyra::createMember(*me->get_x_space());
183 Thyra::assign(zero_.ptr(),0.0);
184 }
185
186 // request only the mass matrix from the physics
187 // Model evaluator builds: alpha*u_dot + beta*F(u) = 0
188 MEB::InArgs<Scalar> inArgs_new = me->createInArgs();
189 inArgs_new.setArgs(inArgs);
190 inArgs_new.set_x_dot(inArgs.get_x_dot());
191 inArgs_new.set_alpha(1.0);
192 inArgs_new.set_beta(0.0);
193
194 // set the one time beta to ensure dirichlet conditions
195 // are correctly included in the mass matrix: do it for
196 // both epetra and Tpetra.
197 if(panzerModel_!=Teuchos::null)
198 panzerModel_->setOneTimeDirichletBeta(1.0);
199#ifdef PANZER_HAVE_EPETRA_STACK
200 else if(panzerEpetraModel_!=Teuchos::null)
201 panzerEpetraModel_->setOneTimeDirichletBeta(1.0);
202#endif
203 else {
204 // assuming the underlying model is a delegator, walk through
205 // the decerator hierarchy until you find a panzer::ME or panzer::EpetraME.
206 // If you don't find one, then throw because you are in a load of trouble anyway!
207 setOneTimeDirichletBeta(1.0,*this->getUnderlyingModel());
208 }
209
210 // set only the mass matrix
211 MEB::OutArgs<Scalar> outArgs = me->createOutArgs();
212 outArgs.set_W_op(mass_);
213
214 // this will fill the mass matrix operator
215 me->evalModel(inArgs_new,outArgs);
216
217 // Teuchos::RCP<const Epetra_CrsMatrix> crsMat = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(Thyra::get_Epetra_Operator(*mass));
218 // EpetraExt::RowMatrixToMatrixMarketFile("expmat.mm",*crsMat);
219
220 if(!massLumping_) {
221 invMassMatrix_ = Thyra::inverse<Scalar>(*me->get_W_factory(),mass_);
222 }
223 else {
224 // build lumped mass matrix (assumes all positive mass entries, does a simple sum)
225 Teuchos::RCP<Thyra::VectorBase<Scalar> > ones = Thyra::createMember(*mass_->domain());
226 Thyra::assign(ones.ptr(),1.0);
227
228 RCP<Thyra::VectorBase<Scalar> > invLumpMass = Thyra::createMember(*mass_->range());
229 Thyra::apply(*mass_,Thyra::NOTRANS,*ones,invLumpMass.ptr());
230 Thyra::reciprocal(*invLumpMass,invLumpMass.ptr());
231
232 invMassMatrix_ = Thyra::diagonal(invLumpMass);
233 }
234}
235
236template<typename Scalar>
239{
240 typedef Thyra::ModelEvaluatorBase MEB;
241
242 MEB::InArgsSetup<Scalar> inArgs(this->getUnderlyingModel()->createInArgs());
243 inArgs.setModelEvalDescription(this->description());
244 inArgs.setSupports(MEB::IN_ARG_alpha,true);
245 inArgs.setSupports(MEB::IN_ARG_beta,true);
246 inArgs.setSupports(MEB::IN_ARG_x_dot,true);
247 prototypeInArgs_ = inArgs;
248
249 MEB::OutArgsSetup<Scalar> outArgs(this->getUnderlyingModel()->createOutArgs());
250 outArgs.setModelEvalDescription(this->description());
251 outArgs.setSupports(MEB::OUT_ARG_W,false);
252 outArgs.setSupports(MEB::OUT_ARG_W_op,false);
253 prototypeOutArgs_ = outArgs;
254}
255
256template<typename Scalar>
258setOneTimeDirichletBeta(double beta,const Thyra::ModelEvaluator<Scalar> & me) const
259{
260 using Teuchos::Ptr;
261 using Teuchos::ptrFromRef;
262 using Teuchos::ptr_dynamic_cast;
263
264 // try to extract base classes that support setOneTimeDirichletBeta
265 Ptr<const panzer::ModelEvaluator<Scalar> > panzerModel = ptr_dynamic_cast<const panzer::ModelEvaluator<Scalar> >(ptrFromRef(me));
266 if(panzerModel!=Teuchos::null) {
267 panzerModel->setOneTimeDirichletBeta(beta);
268 return;
269 }
270 else {
271#ifdef PANZER_HAVE_EPETRA_STACK
272 Ptr<const Thyra::EpetraModelEvaluator> epModel = ptr_dynamic_cast<const Thyra::EpetraModelEvaluator>(ptrFromRef(me));
273 if(epModel!=Teuchos::null) {
274 Ptr<const panzer::ModelEvaluator_Epetra> panzerEpetraModel
275 = ptr_dynamic_cast<const panzer::ModelEvaluator_Epetra>(epModel->getEpetraModel().ptr());
276
277 if(panzerEpetraModel!=Teuchos::null) {
278 panzerEpetraModel->setOneTimeDirichletBeta(beta);
279 return;
280 }
281 }
282#endif
283 }
284
285 // if you get here then the ME is not a panzer::ME or panzer::EpetraME, check
286 // to see if its a delegator
287
288 Ptr<const Thyra::ModelEvaluatorDelegatorBase<Scalar> > delegator
289 = ptr_dynamic_cast<const Thyra::ModelEvaluatorDelegatorBase<Scalar> >(ptrFromRef(me));
290 if(delegator!=Teuchos::null) {
291 setOneTimeDirichletBeta(beta,*delegator->getUnderlyingModel());
292 return;
293 }
294
295 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
296 "panzer::ExplicitModelEvaluator::setOneTimeDirichletBeta can't find a panzer::ME or a panzer::EpetraME. "
297 "The deepest model is also not a delegator. Thus the recursion failed and an exception was generated.");
298}
299
300} // end namespace panzer
301
302#endif
Teuchos::RCP< const panzer::ModelEvaluator< Scalar > > panzerModel_
Access to the panzer model evaluator pointer (thyra version)
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Evaluation of model, applies the mass matrix to the RHS.
void setOneTimeDirichletBeta(double beta, const Thyra::ModelEvaluator< Scalar > &me) const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Build the nominal values, modifies the underlying models in args slightly.
void buildInverseMassMatrix(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs) const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgs() const
Build the out args, modifies the underlying models in args slightly.
Teuchos::RCP< panzer::ModelEvaluator< Scalar > > getPanzerUnderlyingModel()
Get the underlying panzer::ModelEvaluator.
void buildArgsPrototypes()
Build prototype in/out args.
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Build the in args, modifies the underlying models in args slightly.
bool applyMassInverse_
Apply mass matrix inverse within the evaluator.