EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
GLpApp_AdvDiffReactOptModel.cpp
Go to the documentation of this file.
1/*
2//@HEADER
3// ***********************************************************************
4//
5// EpetraExt: Epetra Extended - Linear Algebra Services Package
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 Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//@HEADER
42*/
43
45#include "Epetra_SerialComm.h"
46#include "Epetra_CrsMatrix.h"
47#include "Teuchos_ScalarTraits.hpp"
48#include "Teuchos_VerboseObject.hpp"
49#include "Teuchos_TimeMonitor.hpp"
50
51// 2007/06/14: rabartl: The dependence on Thyra is non-optional right now
52// since I use it to perform a mat-vec that I could not figure out how to do
53// with just epetra. If this becomes a problem then we have to change this
54// later! Just grep for 'Thyra' outside of the ifdef for
55// HAVE_THYRA_EPETRAEXT.
56#include "Thyra_EpetraThyraWrappers.hpp"
57#include "Thyra_VectorStdOps.hpp"
58
59#ifdef HAVE_THYRA_EPETRAEXT
60// For orthogonalization of the basis B_bar
61#include "sillyModifiedGramSchmidt.hpp" // This is just an example!
62#include "Thyra_MultiVectorStdOps.hpp"
63#include "Thyra_SpmdVectorSpaceBase.hpp"
64#endif // HAVE_THYRA_EPETRAEXT
65
66//#define GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
67
68namespace {
69
70Teuchos::RCP<Teuchos::Time>
71 initalizationTimer = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Init"),
72 evalTimer = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval"),
73 p_bar_Timer = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:p_bar"),
74 R_p_bar_Timer = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:R_p_bar"),
75 f_Timer = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:f"),
76 g_Timer = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:g"),
77 W_Timer = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:W"),
78 DfDp_Timer = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:DfDp"),
79 DgDx_Timer = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:DgDx"),
80 DgDp_Timer = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:DgDp");
81
82inline double sqr( const double& s ) { return s*s; }
83
84inline double dot( const Epetra_Vector &x, const Epetra_Vector &y )
85{
86 double dot[1];
87 x.Dot(y,dot);
88 return dot[0];
89}
90
91} // namespace
92
93namespace GLpApp {
94
96 const Teuchos::RCP<const Epetra_Comm> &comm
97 ,const double beta
98 ,const double len_x // Ignored if meshFile is *not* empty
99 ,const double len_y // Ignored if meshFile is *not* empty
100 ,const int local_nx // Ignored if meshFile is *not* empty
101 ,const int local_ny // Ignored if meshFile is *not* empty
102 ,const char meshFile[]
103 ,const int np
104 ,const double x0
105 ,const double p0
106 ,const double reactionRate
107 ,const bool normalizeBasis
108 ,const bool supportDerivatives
109 )
110 : supportDerivatives_(supportDerivatives), np_(np)
111{
112 Teuchos::TimeMonitor initalizationTimerMonitor(*initalizationTimer);
113#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
114 Teuchos::RCP<Teuchos::FancyOStream>
115 out = Teuchos::VerboseObjectBase::getDefaultOStream();
116 Teuchos::OSTab tab(out);
117 *out << "\nEntering AdvDiffReactOptModel::AdvDiffReactOptModel(...) ...\n";
118#endif
119 //
120 // Initalize the data pool object
121 //
122 dat_ = Teuchos::rcp(new GLpApp::GLpYUEpetraDataPool(comm,beta,len_x,len_y,local_nx,local_ny,meshFile,false));
123 //
124 // Get the maps
125 //
126 map_x_ = Teuchos::rcp(new Epetra_Map(dat_->getA()->OperatorDomainMap()));
127 map_p_.resize(Np_);
128 map_p_[p_rx_idx] = Teuchos::rcp(new Epetra_Map(1,1,0,*comm));
129 map_p_bar_ = Teuchos::rcp(new Epetra_Map(dat_->getB()->OperatorDomainMap()));
130 map_f_ = Teuchos::rcp(new Epetra_Map(dat_->getA()->OperatorRangeMap()));
131 map_g_ = Teuchos::rcp(new Epetra_Map(1,1,0,*comm));
132 //
133 // Initialize the basis matrix for p such that p_bar = B_bar * p
134 //
135 if( np_ > 0 ) {
136 //
137 // Create a sine series basis for y (the odd part of a Fourier series basis)
138 //
139 const Epetra_SerialDenseMatrix &ipcoords = *dat_->getipcoords();
140 const Epetra_IntSerialDenseVector &pindx = *dat_->getpindx();
141 Epetra_Map overlapmap(-1,pindx.M(),const_cast<int*>(pindx.A()),1,*comm);
142 const double pi = 2.0 * std::asin(1.0);
143#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
144 *out << "\npi="<<pi<<"\n";
145#endif
146 map_p_[p_bndy_idx] = Teuchos::rcp(new Epetra_Map(np_,np_,0,*comm));
147 B_bar_ = Teuchos::rcp(new Epetra_MultiVector(*map_p_bar_,np_));
148 (*B_bar_)(0)->PutScalar(1.0); // First column is all ones!
149 if( np_ > 1 ) {
150 const int numBndyNodes = map_p_bar_->NumMyElements();
151 const int *bndyNodeGlobalIDs = map_p_bar_->MyGlobalElements();
152 for( int i = 0; i < numBndyNodes; ++i ) {
153#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
154 const int global_id = bndyNodeGlobalIDs[i];
155#endif
156 const int local_id = overlapmap.LID(bndyNodeGlobalIDs[i]);
157 const double x = ipcoords(local_id,0), y = ipcoords(local_id,1);
158 double z = -1.0, L = -1.0;
159 if( x < 1e-10 || len_x - 1e-10 < x ) {
160 z = y;
161 L = len_y;
162 }
163 else {
164 z = x;
165 L = len_x;
166 }
167#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
168 *out << "\ni="<<i<<",global_id="<<global_id<<",local_id="<<local_id<<",x="<<x<<",y="<<y<<",z="<<z<<"\n";
169#endif
170 for( int j = 1; j < np_; ++j ) {
171 (*B_bar_)[j][i] = std::sin(j*pi*z/L);
172#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
173 *out << "\nB("<<i<<","<<j<<")="<<(*B_bar_)[j][i]<<"\n";
174#endif
175 }
176 }
177 if(normalizeBasis) {
178#ifdef HAVE_THYRA_EPETRAEXT
179 //
180 // Use modified Gram-Schmidt to create an orthonormal version of B_bar!
181 //
182 Teuchos::RCP<Thyra::MultiVectorBase<double> >
183 thyra_B_bar = Thyra::create_MultiVector(
184 B_bar_
185 ,Thyra::create_VectorSpace(Teuchos::rcp(new Epetra_Map(*map_p_bar_)))
186 ,Thyra::create_VectorSpace(Teuchos::rcp(new Epetra_Map(*map_p_[p_bndy_idx])))
187 ),
188 thyra_fact_R;
189 Thyra::sillyModifiedGramSchmidt(thyra_B_bar.ptr(), Teuchos::outArg(thyra_fact_R));
190 Thyra::scale(double(numBndyNodes)/double(np_), thyra_B_bar.ptr()); // Each row should sum to around one!
191 // We just discard the "R" factory thyra_fact_R ...
192#else // HAVE_THYRA_EPETRAEXT
193 TEUCHOS_TEST_FOR_EXCEPTION(
194 true,std::logic_error
195 ,"Error, can not normalize basis since we do not have Thyra support enabled!"
196 );
197#endif // HAVE_EPETRAEXT_THYRA
198 }
199 }
200#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
201 *out << "\nB_bar =\n\n";
202 B_bar_->Print(Teuchos::OSTab(out).o());
203#endif
204 }
205 else {
206 // B_bar = I implicitly!
208 }
209 //
210 // Create vectors
211 //
212 x0_ = Teuchos::rcp(new Epetra_Vector(*map_x_));
213 //xL_ = Teuchos::rcp(new Epetra_Vector(*map_x_));
214 //xU_ = Teuchos::rcp(new Epetra_Vector(*map_x_));
215 p0_.resize(Np_);
216 p0_[p_bndy_idx] = Teuchos::rcp(new Epetra_Vector(*map_p_[p_bndy_idx]));
217 p0_[p_rx_idx] = Teuchos::rcp(new Epetra_Vector(*map_p_[p_rx_idx]));
218 pL_.resize(Np_);
219 pU_.resize(Np_);
220 //pL_ = Teuchos::rcp(new Epetra_Vector(*map_p_));
221 //pU_ = Teuchos::rcp(new Epetra_Vector(*map_p_));
222 //
223 // Initialize the vectors
224 //
225 x0_->PutScalar(x0);
226 p0_[p_bndy_idx]->PutScalar(p0);
227 p0_[p_rx_idx]->PutScalar(reactionRate);
228 //
229 // Initialize the graph for W
230 //
231 dat_->computeNpy(x0_);
232 //if(dumpAll) { *out << "\nA =\n"; { Teuchos::OSTab tab(out); dat_->getA()->Print(*out); } }
233 //if(dumpAll) { *out << "\nNpy =\n"; { Teuchos::OSTab tab(out); dat_->getNpy()->Print(*out); } }
234 W_graph_ = Teuchos::rcp(new Epetra_CrsGraph(dat_->getA()->Graph())); // Assume A and Npy have same graph!
235 //
236 // Get default objective matching vector q
237 //
238 q_ = Teuchos::rcp(new Epetra_Vector(*(*dat_->getq())(0))); // From Epetra_FEVector to Epetra_Vector!
239 //
240 isInitialized_ = true;
241#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
242 *out << "\nLeaving AdvDiffReactOptModel::AdvDiffReactOptModel(...) ...\n";
243#endif
244}
245
246void AdvDiffReactOptModel::set_q( Teuchos::RCP<const Epetra_Vector> const& q )
247{
248 q_ = q;
249#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
250 Teuchos::RCP<Teuchos::FancyOStream>
251 dout = Teuchos::VerboseObjectBase::getDefaultOStream();
252 Teuchos::OSTab dtab(dout);
253 *dout << "\nIn AdvDiffReactOptModel::set_q(q): q =\n\n";
254 q_->Print(Teuchos::OSTab(dout).o());
255#endif
256}
257
258Teuchos::RCP<GLpApp::GLpYUEpetraDataPool>
260{
261 return dat_;
262}
263
264Teuchos::RCP<const Epetra_MultiVector>
266{
267 return B_bar_;
268}
269
270// Overridden from EpetraExt::ModelEvaluator
271
272Teuchos::RCP<const Epetra_Map>
274{
275 return map_x_;
276}
277
278Teuchos::RCP<const Epetra_Map>
280{
281 return map_x_;
282}
283
284Teuchos::RCP<const Epetra_Map>
286{
287 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=Np_));
288 return map_p_[l];
289}
290
291Teuchos::RCP<const Epetra_Map>
293{
294 TEUCHOS_TEST_FOR_EXCEPT(j!=0);
295 return map_g_;
296}
297
298Teuchos::RCP<const Epetra_Vector>
300{
301 return x0_;
302}
303
304Teuchos::RCP<const Epetra_Vector>
306{
307 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=Np_));
308 return p0_[l];
309}
310
311Teuchos::RCP<const Epetra_Vector>
313{
314 return xL_;
315}
316
317Teuchos::RCP<const Epetra_Vector>
319{
320 return xU_;
321}
322
323Teuchos::RCP<const Epetra_Vector>
325{
326 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=Np_));
327 return pL_[l];
328}
329
330Teuchos::RCP<const Epetra_Vector>
332{
333 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=Np_));
334 return pU_[l];
335}
336
337Teuchos::RCP<Epetra_Operator>
339{
340 return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
341}
342
343Teuchos::RCP<Epetra_Operator>
345{
346 TEUCHOS_TEST_FOR_EXCEPT(l!=0);
347 return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,dat_->getB()->Graph()));
348 // See DfDp in evalModel(...) below for details
349}
350
353{
354 InArgsSetup inArgs;
355 inArgs.setModelEvalDescription(this->description());
356 inArgs.set_Np(2);
357 inArgs.setSupports(IN_ARG_x,true);
358 return inArgs;
359}
360
363{
364 OutArgsSetup outArgs;
365 outArgs.setModelEvalDescription(this->description());
366 outArgs.set_Np_Ng(2,1);
367 outArgs.setSupports(OUT_ARG_f,true);
368 outArgs.setSupports(OUT_ARG_W,true);
369 outArgs.set_W_properties(
373 ,true // supportsAdjoint
374 )
375 );
377 outArgs.setSupports(
379 ,( np_ > 0
382 )
383 );
384 outArgs.set_DfDp_properties(
388 ,true // supportsAdjoint
389 )
390 );
392 outArgs.set_DgDx_properties(
396 ,true // supportsAdjoint
397 )
398 );
400 outArgs.set_DgDp_properties(
404 ,true // supportsAdjoint
405 )
406 );
407 }
408 return outArgs;
409}
410
411void AdvDiffReactOptModel::evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const
412{
413 Teuchos::TimeMonitor evalTimerMonitor(*evalTimer);
414#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
415 Teuchos::RCP<Teuchos::FancyOStream>
416 dout = Teuchos::VerboseObjectBase::getDefaultOStream();
417 Teuchos::OSTab dtab(dout);
418 *dout << "\nEntering AdvDiffReactOptModel::evalModel(...) ...\n";
419#endif
420 using Teuchos::dyn_cast;
421 using Teuchos::rcp_dynamic_cast;
422 //
423 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
424 const bool trace = ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_LOW) );
425 const bool dumpAll = ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_EXTREME) );
426 //
427 Teuchos::OSTab tab(out);
428 if(out.get() && trace) *out << "\n*** Entering AdvDiffReactOptModel::evalModel(...) ...\n";
429 //
430 // Get the input arguments
431 //
432 const Epetra_Vector *p_in = inArgs.get_p(p_bndy_idx).get();
433 const Epetra_Vector &p = (p_in ? *p_in : *p0_[p_bndy_idx]);
434 const Epetra_Vector *rx_vec_in = inArgs.get_p(p_rx_idx).get();
435 const double reactionRate = ( rx_vec_in ? (*rx_vec_in)[0] : (*p0_[p_rx_idx])[0] );
436 const Epetra_Vector &x = *inArgs.get_x();
437#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
438 *dout << "\nx =\n\n";
439 x.Print(Teuchos::OSTab(dout).o());
440 *dout << "\np =\n\n";
441 p.Print(Teuchos::OSTab(dout).o());
442#endif
443 //
444 // Get the output arguments
445 //
446 Epetra_Vector *f_out = outArgs.get_f().get();
447 Epetra_Vector *g_out = outArgs.get_g(0).get();
448 Epetra_Operator *W_out = outArgs.get_W().get();
449 Derivative DfDp_out;
450 Epetra_MultiVector *DgDx_trans_out = 0;
451 Epetra_MultiVector *DgDp_trans_out = 0;
453 DfDp_out = outArgs.get_DfDp(0);
454 DgDx_trans_out = get_DgDx_mv(0,outArgs,DERIV_TRANS_MV_BY_ROW).get();
455 DgDp_trans_out = get_DgDp_mv(0,0,outArgs,DERIV_TRANS_MV_BY_ROW).get();
456 }
457 //
458 // Precompute some shared quantities
459 //
460 // p_bar = B_bar*p
461 Teuchos::RCP<const Epetra_Vector> p_bar;
462 if(np_ > 0) {
463 Teuchos::TimeMonitor p_bar_TimerMonitor(*p_bar_Timer);
464 Teuchos::RCP<Epetra_Vector> _p_bar;
465 _p_bar = Teuchos::rcp(new Epetra_Vector(*map_p_bar_));
466 _p_bar->Multiply('N','N',1.0,*B_bar_,p,0.0);
467 p_bar = _p_bar;
468 }
469 else {
470 p_bar = Teuchos::rcp(&p,false);
471 }
472 // R_p_bar = R*p_bar = R*(B_bar*p)
473 Teuchos::RCP<const Epetra_Vector> R_p_bar;
474 if( g_out || DgDp_trans_out ) {
475 Teuchos::TimeMonitor R_p_bar_TimerMonitor(*R_p_bar_Timer);
476 Teuchos::RCP<Epetra_Vector>
477 _R_p_bar = Teuchos::rcp(new Epetra_Vector(*map_p_bar_));
478 dat_->getR()->Multiply(false,*p_bar,*_R_p_bar);
479 R_p_bar = _R_p_bar;
480 }
481 //
482 // Compute the functions
483 //
484 if(f_out) {
485 //
486 // f = A*x + reationRate*Ny(x) + B*(B_bar*p)
487 //
488 Teuchos::TimeMonitor f_TimerMonitor(*f_Timer);
489 Epetra_Vector &f = *f_out;
491 dat_->getA()->Multiply(false,x,Ax);
492 f.Update(1.0,Ax,0.0);
493 if(reactionRate!=0.0) {
494 dat_->computeNy(Teuchos::rcp(&x,false));
495 f.Update(reactionRate,*dat_->getNy(),1.0);
496 }
498 dat_->getB()->Multiply(false,*p_bar,Bp);
499 f.Update(1.0,Bp,-1.0,*dat_->getb(),1.0);
500 }
501 if(g_out) {
502 //
503 // g = 0.5 * (x-q)'*H*(x-q) + 0.5*regBeta*(B_bar*p)'*R*(B_bar*p)
504 //
505 Teuchos::TimeMonitor g_TimerMonitor(*g_Timer);
506 Epetra_Vector &g = *g_out;
507 Epetra_Vector xq(x);
508 xq.Update(-1.0, *q_, 1.0);
509 Epetra_Vector Hxq(x);
510 dat_->getH()->Multiply(false,xq,Hxq);
511 g[0] = 0.5*dot(xq,Hxq) + 0.5*dat_->getbeta()*dot(*p_bar,*R_p_bar);
512#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
513 *dout << "q =\n\n";
514 q_->Print(Teuchos::OSTab(dout).o());
515 *dout << "x =\n\n";
516 x.Print(Teuchos::OSTab(dout).o());
517 *dout << "x-q =\n\n";
518 xq.Print(Teuchos::OSTab(dout).o());
519 *dout << "H*(x-q) =\n\n";
520 Hxq.Print(Teuchos::OSTab(dout).o());
521 *dout << "0.5*dot(xq,Hxq) = " << (0.5*dot(xq,Hxq)) << "\n";
522 *dout << "0.5*regBeta*(B_bar*p)'*R*(B_bar*p) = " << (0.5*dat_->getbeta()*dot(*p_bar,*R_p_bar)) << "\n";
523#endif
524 }
525 if(W_out) {
526 //
527 // W = A + reationRate*Npy(x)
528 //
529 Teuchos::TimeMonitor W_TimerMonitor(*W_Timer);
530 Epetra_CrsMatrix &DfDx = dyn_cast<Epetra_CrsMatrix>(*W_out);
531 if(reactionRate!=0.0)
532 dat_->computeNpy(Teuchos::rcp(&x,false));
533 Teuchos::RCP<Epetra_CrsMatrix>
534 dat_A = dat_->getA(),
535 dat_Npy = dat_->getNpy();
536 const int numMyRows = dat_A->NumMyRows();
537 for( int i = 0; i < numMyRows; ++i ) {
538 int dat_A_num_row_entries=0; double *dat_A_row_vals=0; int *dat_A_row_inds=0;
539 dat_A->ExtractMyRowView(i,dat_A_num_row_entries,dat_A_row_vals,dat_A_row_inds);
540 int DfDx_num_row_entries=0; double *DfDx_row_vals=0; int *DfDx_row_inds=0;
541 DfDx.ExtractMyRowView(i,DfDx_num_row_entries,DfDx_row_vals,DfDx_row_inds);
542#ifdef TEUCHOS_DEBUG
543 TEUCHOS_TEST_FOR_EXCEPT(DfDx_num_row_entries!=dat_A_num_row_entries);
544#endif
545 if(reactionRate!=0.0) {
546 int dat_Npy_num_row_entries=0; double *dat_Npy_row_vals=0; int *dat_Npy_row_inds=0;
547 dat_Npy->ExtractMyRowView(i,dat_Npy_num_row_entries,dat_Npy_row_vals,dat_Npy_row_inds);
548#ifdef TEUCHOS_DEBUG
549 TEUCHOS_TEST_FOR_EXCEPT(dat_A_num_row_entries!=dat_Npy_num_row_entries);
550#endif
551 for(int k = 0; k < DfDx_num_row_entries; ++k) {
552#ifdef TEUCHOS_DEBUG
553 TEUCHOS_TEST_FOR_EXCEPT(dat_A_row_inds[k]!=dat_Npy_row_inds[k]||dat_A_row_inds[k]!=DfDx_row_inds[k]);
554#endif
555 DfDx_row_vals[k] = dat_A_row_vals[k] + reactionRate * dat_Npy_row_vals[k];
556 }
557 }
558 else {
559 for(int k = 0; k < DfDx_num_row_entries; ++k) {
560#ifdef TEUCHOS_DEBUG
561 TEUCHOS_TEST_FOR_EXCEPT(dat_A_row_inds[k]!=DfDx_row_inds[k]);
562#endif
563 DfDx_row_vals[k] = dat_A_row_vals[k];
564 }
565 }
566 }
567 }
568 if(!DfDp_out.isEmpty()) {
569 if(out.get() && trace) *out << "\nComputing DfDp ...\n";
570 //
571 // DfDp = B*B_bar
572 //
573 Teuchos::TimeMonitor DfDp_TimerMonitor(*DfDp_Timer);
574 Epetra_CrsMatrix *DfDp_op = NULL;
575 Epetra_MultiVector *DfDp_mv = NULL;
576 if(out.get() && dumpAll)
577 { *out << "\nB =\n"; { Teuchos::OSTab tab(out); dat_->getB()->Print(*out); } }
578 if(DfDp_out.getLinearOp().get()) {
579 DfDp_op = &dyn_cast<Epetra_CrsMatrix>(*DfDp_out.getLinearOp());
580 }
581 else {
582 DfDp_mv = &*DfDp_out.getDerivativeMultiVector().getMultiVector();
583 DfDp_mv->PutScalar(0.0);
584 }
585 Teuchos::RCP<Epetra_CrsMatrix>
586 dat_B = dat_->getB();
587 if(np_ > 0) {
588 //
589 // We only support a Multi-vector form when we have a non-idenity basis
590 // matrix B_bar for p!
591 //
592 TEUCHOS_TEST_FOR_EXCEPT(DfDp_mv==NULL);
593 dat_B->Multiply(false,*B_bar_,*DfDp_mv);
594 }
595 else {
596 //
597 // Note: We copy from B every time in order to be safe. Note that since
598 // the client knows that B is constant (sense we told them so in
599 // createOutArgs()) then it should only compute this matrix once and keep
600 // it if it is smart.
601 //
602 // Note: We support both the CrsMatrix and MultiVector form of this object
603 // to make things easier for the client.
604 //
605 if(DfDp_op) {
606 const int numMyRows = dat_B->NumMyRows();
607 for( int i = 0; i < numMyRows; ++i ) {
608 int dat_B_num_row_entries=0; double *dat_B_row_vals=0; int *dat_B_row_inds=0;
609 dat_B->ExtractMyRowView(i,dat_B_num_row_entries,dat_B_row_vals,dat_B_row_inds);
610 int DfDp_num_row_entries=0; double *DfDp_row_vals=0; int *DfDp_row_inds=0;
611 DfDp_op->ExtractMyRowView(i,DfDp_num_row_entries,DfDp_row_vals,DfDp_row_inds);
612#ifdef TEUCHOS_DEBUG
613 TEUCHOS_TEST_FOR_EXCEPT(DfDp_num_row_entries!=dat_B_num_row_entries);
614#endif
615 for(int k = 0; k < DfDp_num_row_entries; ++k) {
616#ifdef TEUCHOS_DEBUG
617 TEUCHOS_TEST_FOR_EXCEPT(dat_B_row_inds[k]!=DfDp_row_inds[k]);
618#endif
619 DfDp_row_vals[k] = dat_B_row_vals[k];
620 }
621 // ToDo: The above code should be put in a utility function called copyValues(...)!
622 }
623 }
624 else if(DfDp_mv) {
625 // We must do a mat-vec to get this since we can't just copy out the
626 // matrix entries since the domain map may be different from the
627 // column map! I learned this the very very hard way! I am using
628 // Thyra wrappers here since I can't figure out for the life of me how
629 // to do this cleanly with Epetra alone!
630 Teuchos::RCP<Epetra_Vector>
631 etaVec = Teuchos::rcp(new Epetra_Vector(*map_p_bar_));
632 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
633 space_p_bar = Thyra::create_VectorSpace(Teuchos::rcp(new Epetra_Map(*map_p_bar_)));
634 Teuchos::RCP<Thyra::VectorBase<double> >
635 thyra_etaVec = Thyra::create_Vector(etaVec,space_p_bar);
636 for( int i = 0; i < map_p_bar_->NumGlobalElements(); ++i ) {
637 Thyra::assign(thyra_etaVec.ptr(), 0.0);
638 Thyra::set_ele(i, 1.0, thyra_etaVec.ptr());
639 dat_B->Multiply(false, *etaVec, *(*DfDp_mv)(i));
640 };
641 }
642 }
643#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
644 if(DfDp_op) {
645 *dout << "\nDfDp_op =\n\n";
646 DfDp_op->Print(Teuchos::OSTab(dout).o());
647 }
648 if(DfDp_mv) {
649 *dout << "\nDfDp_mv =\n\n";
650 DfDp_mv->Print(Teuchos::OSTab(dout).o());
651 }
652#endif
653 }
654 if(DgDx_trans_out) {
655 //
656 // DgDx' = H*(x-q)
657 //
658 Teuchos::TimeMonitor DgDx_TimerMonitor(*DgDx_Timer);
659 Epetra_Vector &DgDx_trans = *(*DgDx_trans_out)(0);
660 Epetra_Vector xq(x);
661 xq.Update(-1.0,*q_,1.0);
662 dat_->getH()->Multiply(false,xq,DgDx_trans);
663#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
664 *dout << "q =\n\n";
665 q_->Print(Teuchos::OSTab(dout).o());
666 *dout << "x =\n\n";
667 x.Print(Teuchos::OSTab(dout).o());
668 *dout << "x-q =\n\n";
669 xq.Print(Teuchos::OSTab(dout).o());
670 *dout << "DgDx_trans = H*(x-q) =\n\n";
671 DgDx_trans.Print(Teuchos::OSTab(dout).o());
672#endif
673 }
674 if(DgDp_trans_out) {
675 //
676 // DgDp' = regBeta*B_bar'*(R*(B_bar*p))
677 //
678 Teuchos::TimeMonitor DgDp_TimerMonitor(*DgDp_Timer);
679 Epetra_Vector &DgDp_trans = *(*DgDp_trans_out)(0);
680 if(np_ > 0) {
681 DgDp_trans.Multiply('T','N',dat_->getbeta(),*B_bar_,*R_p_bar,0.0);
682 }
683 else {
684 DgDp_trans.Update(dat_->getbeta(),*R_p_bar,0.0);
685 }
686#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
687 *dout << "\nR_p_bar =\n\n";
688 R_p_bar->Print(Teuchos::OSTab(dout).o());
689 if(B_bar_.get()) {
690 *dout << "\nB_bar =\n\n";
691 B_bar_->Print(Teuchos::OSTab(dout).o());
692 }
693 *dout << "\nDgDp_trans =\n\n";
694 DgDp_trans.Print(Teuchos::OSTab(dout).o());
695#endif
696 }
697 if(out.get() && trace) *out << "\n*** Leaving AdvDiffReactOptModel::evalModel(...) ...\n";
698#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
699 *dout << "\nLeaving AdvDiffReactOptModel::evalModel(...) ...\n";
700#endif
701}
702
703} // namespace GLpApp
Teuchos::RCP< Epetra_MultiVector > getMultiVector() const
Simple aggregate class that stores a derivative object as a general linear operator or as a multi-vec...
DerivativeMultiVector getDerivativeMultiVector() const
Teuchos::RCP< Epetra_Operator > getLinearOp() const
void setSupports(EInArgsMembers arg, bool supports=true)
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< const Epetra_Vector > get_p(int l) const
Teuchos::RCP< const Epetra_Vector > get_x() const
Set solution vector Taylor polynomial.
void setModelEvalDescription(const std::string &modelEvalDescription)
void set_W_properties(const DerivativeProperties &properties)
void set_DfDp_properties(int l, const DerivativeProperties &properties)
void set_DgDp_properties(int j, int l, const DerivativeProperties &properties)
void set_DgDx_properties(int j, const DerivativeProperties &properties)
void setSupports(EOutArgsMembers arg, bool supports=true)
Teuchos::RCP< Epetra_Operator > get_W() const
Evaluation< Epetra_Vector > get_g(int j) const
Get g(j) where 0 <= j && j < this->Ng().
Evaluation< Epetra_Vector > get_f() const
int LID(int GID) const
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
virtual void Print(std::ostream &os) const
const int * A() const
int Dot(const Epetra_MultiVector &A, double *Result) const
virtual void Print(std::ostream &os) const
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int PutScalar(double ScalarConstant)
Teuchos::RCP< Epetra_Operator > create_W() const
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
\breif .
Teuchos::RCP< Epetra_MultiVector > B_bar_
Teuchos::RCP< const Epetra_Map > map_x_
Teuchos::RCP< const Epetra_Map > map_p_bar_
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Teuchos::RCP< const Epetra_Vector > get_p_upper_bounds(int l) const
Teuchos::RCP< GLpApp::GLpYUEpetraDataPool > dat_
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
\breif .
Teuchos::RCP< const Epetra_Vector > q_
Teuchos::RCP< const Epetra_Map > map_f_
Teuchos::RCP< Epetra_CrsGraph > W_graph_
Teuchos::RCP< GLpApp::GLpYUEpetraDataPool > getDataPool()
Teuchos::RCP< const Epetra_Vector > get_x_upper_bounds() const
Teuchos::RCP< const Epetra_Vector > get_p_lower_bounds(int l) const
Teuchos::RCP< Epetra_Operator > create_DfDp_op(int l) const
Teuchos::RCP< const Epetra_MultiVector > get_B_bar() const
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Teuchos::RCP< const Epetra_Map > get_f_map() const
void set_q(Teuchos::RCP< const Epetra_Vector > const &q)
Teuchos::RCP< const Epetra_Map > map_g_
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
AdvDiffReactOptModel(const Teuchos::RCP< const Epetra_Comm > &comm, const double beta, const double len_x, const double len_y, const int local_nx, const int local_ny, const char meshFile[], const int np, const double x0, const double p0, const double reactionRate, const bool normalizeBasis, const bool supportDerivatives)
Constructor.
Teuchos::RCP< const Epetra_Vector > get_x_lower_bounds() const
Teuchos::RCP< const Epetra_Map > get_x_map() const