Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultMultiPeriodModelEvaluator.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) 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 Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef THYRA_DEFAULT_MULTI_PERIOD_MODEL_EVALUATOR_HPP
43#define THYRA_DEFAULT_MULTI_PERIOD_MODEL_EVALUATOR_HPP
44
45
46#include "Thyra_ModelEvaluatorDefaultBase.hpp"
47#include "Thyra_DefaultProductVectorSpace.hpp"
48#include "Thyra_DefaultBlockedTriangularLinearOpWithSolveFactory.hpp" // Default implementation
49#include "Thyra_DefaultBlockedLinearOp.hpp"
50#include "Thyra_ModelEvaluatorDelegatorBase.hpp"
51#include "Thyra_ProductVectorBase.hpp"
52#include "Teuchos_implicit_cast.hpp"
53#include "Teuchos_AbstractFactory.hpp" // Interface
54#include "Teuchos_AbstractFactoryStd.hpp" // Implementation
55
56
57namespace Thyra {
58
59
130template<class Scalar>
132 : virtual public ModelEvaluatorDefaultBase<Scalar>
133{
134public:
135
138
141
144 const int N,
145 const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
146 const Array<int> &z_indexes,
147 const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
148 const int g_index,
149 const Array<Scalar> g_weights,
150 const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space = Teuchos::null,
151 const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space = Teuchos::null,
152 const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory = Teuchos::null
153 );
154
213 void initialize(
214 const int N,
215 const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
216 const Array<int> &z_indexes,
217 const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
218 const int g_index,
219 const Array<Scalar> g_weights,
220 const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space = Teuchos::null,
221 const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space = Teuchos::null,
222 const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory = Teuchos::null
223 );
224
234 void reset_z(
235 const Array<Array<RCP<const VectorBase<Scalar> > > > &z
236 );
237
239
242
244 int Np() const;
246 int Ng() const;
274 void reportFinalPoint(
275 const ModelEvaluatorBase::InArgs<Scalar> &finalPoint,
276 const bool wasSolved
277 );
278
280
281private:
282
283
286
288 RCP<LinearOpBase<Scalar> > create_DfDp_op_impl(int l) const;
290 RCP<LinearOpBase<Scalar> > create_DgDx_dot_op_impl(int j) const;
292 RCP<LinearOpBase<Scalar> > create_DgDx_op_impl(int j) const;
294 RCP<LinearOpBase<Scalar> > create_DgDp_op_impl(int j, int l) const;
296 ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
298 void evalModelImpl(
301 ) const;
302
304
305private:
306
307 // //////////////////////////
308 // Private types
309
312
313 // /////////////////////////
314 // Private data members
315
316 RCP<ModelEvaluator<Scalar> > periodModel_;
317 Array<RCP<ModelEvaluator<Scalar> > > periodModels_;
318 Array<int> z_indexes_;
319 Array<int> period_l_map_;
320 z_t z_; // size == N
321 int g_index_;
322 g_weights_t g_weights_; // size == N
326 int Np_;
327 int Ng_;
331
332 // /////////////////////////
333 // Private member functions
334
335 void set_z_indexes_and_create_period_l_map( const Array<int> &z_indexes );
336
337 void wrapNominalValuesAndBounds();
338
339 static
341 createProductVector(
342 const RCP<const ProductVectorSpaceBase<Scalar> > &prodVecSpc
343 );
344
345 // Return the index of a "free" parameter in the period model given its
346 // index in this mulit-period model.
347 int period_l(const int l) const
348 {
349 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= l && l < Np_) );
350 return period_l_map_[l];
351 }
352
353 int numPeriodZs() const { return z_indexes_.size(); }
354
355 int N() const { return x_bar_space_->numBlocks(); }
356
357};
358
359
360// /////////////////////////////////
361// Implementations
362
363
364// Constructors/Intializers/Accessors
365
366
367template<class Scalar>
369 :g_index_(-1), Np_(-1), Ng_(-1)
370{}
371
372
373template<class Scalar>
375 const int N,
376 const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
377 const Array<int> &z_indexes,
378 const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
379 const int g_index,
380 const Array<Scalar> g_weights,
381 const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space,
382 const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space,
383 const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory
384 )
385 :g_index_(-1), Np_(-1), Ng_(-1)
386{
388 N, periodModels, z_indexes, z, g_index, g_weights,
389 x_bar_space, f_bar_space, W_bar_factory
390 );
391}
392
393
394template<class Scalar>
396 const int N,
397 const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
398 const Array<int> &z_indexes,
399 const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
400 const int g_index,
401 const Array<Scalar> g_weights,
402 const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space,
403 const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space,
404 const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory
405 )
406{
407
410 typedef ModelEvaluatorBase MEB;
411
412 TEUCHOS_TEST_FOR_EXCEPT( N <= 0 );
413 TEUCHOS_TEST_FOR_EXCEPT( implicit_cast<int>(periodModels.size()) != N );
414 MEB::InArgs<Scalar> periodInArgs = periodModels[0]->createInArgs();
415 MEB::OutArgs<Scalar> periodOutArgs = periodModels[0]->createOutArgs();
416 for ( int k = 0; k < implicit_cast<int>(z_indexes.size()); ++k ) {
417 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= z_indexes[k] && z_indexes[k] < periodInArgs.Np() ) );
418 }
419 TEUCHOS_TEST_FOR_EXCEPT( implicit_cast<int>(z.size())!=N );
420 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= g_index && g_index < periodOutArgs.Ng() ) );
421 TEUCHOS_TEST_FOR_EXCEPT( implicit_cast<int>(g_weights.size())!=N );
422 // ToDo: Assert that the different period models are compatible!
423
424 Np_ = periodInArgs.Np() - z_indexes.size();
425
426 Ng_ = 1;
427
428 set_z_indexes_and_create_period_l_map(z_indexes);
429
430 periodModel_ = periodModels[0]; // Assume basic structure!
431
432 periodModels_ = periodModels;
433
434 z_.resize(N);
435 reset_z(z);
436
437 g_index_ = g_index;
438 g_weights_ = g_weights;
439
440 if ( periodInArgs.supports(MEB::IN_ARG_x) ) {
441 if( !is_null(x_bar_space) ) {
442 TEUCHOS_TEST_FOR_EXCEPT(!(x_bar_space->numBlocks()==N));
443 // ToDo: Check the constituent spaces more carefully against models[]->get_x_space().
444 x_bar_space_ = x_bar_space;
445 }
446 else {
447 x_bar_space_ = productVectorSpace(periodModel_->get_x_space(),N);
448 // ToDo: Update to build from different models for different x spaces!
449 }
450 }
451
452 if ( periodOutArgs.supports(MEB::OUT_ARG_f) ) {
453 if( !is_null(f_bar_space) ) {
454 TEUCHOS_TEST_FOR_EXCEPT(!(f_bar_space->numBlocks()==N));
455 // ToDo: Check the constituent spaces more carefully against models[]->get_f_space().
456 f_bar_space_ = f_bar_space;
457 }
458 else {
459 f_bar_space_ = productVectorSpace(periodModel_->get_f_space(),N);
460 // ToDo: Update to build from different models for different f spaces!
461 }
462 }
463
464 if ( periodOutArgs.supports(MEB::OUT_ARG_W) ) {
465 if ( !is_null(W_bar_factory) ) {
466 // ToDo: Check the compatability of the W_bar objects created using this object!
467 W_bar_factory_ = W_bar_factory;
468 }
469 else {
470 W_bar_factory_ =
471 defaultBlockedTriangularLinearOpWithSolveFactory<Scalar>(
472 periodModel_->get_W_factory() );
473 }
474 }
475
476 wrapNominalValuesAndBounds();
477
478}
479
480
481template<class Scalar>
483 const Array<Array<RCP<const VectorBase<Scalar> > > > &z
484 )
485{
486
488
489 const int N = z_.size();
490
491#ifdef TEUCHOS_DEBUG
492 TEUCHOS_TEST_FOR_EXCEPT( N == 0 && "Error, must have called initialize() first!" );
493 TEUCHOS_TEST_FOR_EXCEPT( implicit_cast<int>(z.size()) != N );
494#endif
495
496 for( int i = 0; i < N; ++i ) {
497 const Array<RCP<const VectorBase<Scalar> > > &z_i = z[i];
498#ifdef TEUCHOS_DEBUG
499 TEUCHOS_TEST_FOR_EXCEPT( z_i.size() != z_indexes_.size() );
500#endif
501 z_[i] = z_i;
502 }
503
504}
505
506
507// Public functions overridden from ModelEvaulator
508
509
510template<class Scalar>
512{
513 return Np_;
514}
515
516
517template<class Scalar>
519{
520 return Ng_;
521}
522
523
524template<class Scalar>
527{
528 return x_bar_space_;
529}
530
531
532template<class Scalar>
535{
536 return f_bar_space_;
537}
538
539
540template<class Scalar>
543{
544 return periodModel_->get_p_space(period_l(l));
545}
546
547
548template<class Scalar>
551{
552 return periodModel_->get_p_names(period_l(l));
553}
554
555
556template<class Scalar>
559{
561 return periodModel_->get_g_space(g_index_);
562}
563
564
565template<class Scalar>
568{
569 return periodModel_->get_g_names(j);
570}
571
572
573template<class Scalar>
576{
577 return nominalValues_;
578}
579
580
581template<class Scalar>
584{
585 return lowerBounds_;
586}
587
588
589template<class Scalar>
592{
593 return upperBounds_;
594}
595
596
597template<class Scalar>
600{
601 // Set up the block structure ready to have the blocks filled!
603 W_op_bar = defaultBlockedLinearOp<Scalar>();
604 W_op_bar->beginBlockFill(f_bar_space_,x_bar_space_);
605 const int N = x_bar_space_->numBlocks();
606 for ( int i = 0; i < N; ++i ) {
607 W_op_bar->setNonconstBlock( i, i, periodModel_->create_W_op() );
608 }
609 W_op_bar->endBlockFill();
610 return W_op_bar;
611}
612
613
614template<class Scalar>
617{
618 return Teuchos::null;
619}
620
621
622template<class Scalar>
625{
626 return W_bar_factory_;
627}
628
629
630template<class Scalar>
633{
634 typedef ModelEvaluatorBase MEB;
635 MEB::InArgs<Scalar> periodInArgs = periodModel_->createInArgs();
636 MEB::InArgsSetup<Scalar> inArgs;
637 inArgs.setModelEvalDescription(this->description());
638 inArgs.set_Np(Np_);
639 inArgs.setSupports( MEB::IN_ARG_x, periodInArgs.supports(MEB::IN_ARG_x) );
640 return inArgs;
641}
642
643
644template<class Scalar>
646 const ModelEvaluatorBase::InArgs<Scalar> &finalPoint
647 ,const bool wasSolved
648 )
649{
650 // We are just going to ignore the final point here. It is not clear how to
651 // report a "final" point back to the underlying *periodModel_ object since
652 // we have so many different "points" that we could return (i.e. one for
653 // each period). I guess we could report back the final parameter values
654 // (other than the z parameter) but there are multiple states x[i] and
655 // period parameters z[i] that we can report back.
656}
657
658
659// Public functions overridden from ModelEvaulatorDefaultBase
660
661
662template<class Scalar>
665{
666 TEUCHOS_TEST_FOR_EXCEPT("This class does not support DfDp(l) as a linear operator yet.");
667 return Teuchos::null;
668}
669
670
671template<class Scalar>
672RCP<LinearOpBase<Scalar> >
673DefaultMultiPeriodModelEvaluator<Scalar>::create_DgDx_dot_op_impl(int j) const
674{
675 TEUCHOS_TEST_FOR_EXCEPT("This class does not support DgDx_dot(j) as a linear operator yet.");
676 return Teuchos::null;
677}
678
679
680template<class Scalar>
681RCP<LinearOpBase<Scalar> >
682DefaultMultiPeriodModelEvaluator<Scalar>::create_DgDx_op_impl(int j) const
683{
684 TEUCHOS_TEST_FOR_EXCEPT("This class does not support DgDx(j) as a linear operator yet.");
685 return Teuchos::null;
686}
687
688
689template<class Scalar>
690RCP<LinearOpBase<Scalar> >
691DefaultMultiPeriodModelEvaluator<Scalar>::create_DgDp_op_impl(int j, int l) const
692{
693 TEUCHOS_TEST_FOR_EXCEPT("This class does not support DgDp(j,l) as a linear operator yet.");
694 return Teuchos::null;
695}
696
697
698template<class Scalar>
699ModelEvaluatorBase::OutArgs<Scalar>
700DefaultMultiPeriodModelEvaluator<Scalar>::createOutArgsImpl() const
701{
702
703 typedef ModelEvaluatorBase MEB;
704
705 MEB::OutArgs<Scalar> periodOutArgs = periodModel_->createOutArgs();
706 MEB::OutArgsSetup<Scalar> outArgs;
707
708 outArgs.setModelEvalDescription(this->description());
709
710 outArgs.set_Np_Ng(Np_,Ng_);
711
712 // f
713 if (periodOutArgs.supports(MEB::OUT_ARG_f) ) {
714 outArgs.setSupports(MEB::OUT_ARG_f);
715 }
716
717 // W_op
718 if (periodOutArgs.supports(MEB::OUT_ARG_W_op) ) {
719 outArgs.setSupports(MEB::OUT_ARG_W_op);
720 outArgs.set_W_properties(periodOutArgs.get_W_properties());
721 }
722 // Note: We will not directly support the LOWSB form W as we will let the
723 // default base class handle this given our W_factory!
724
725 // DfDp(l)
726 for ( int l = 0; l < Np_; ++l ) {
727 const int period_l = this->period_l(l);
728 const MEB::DerivativeSupport period_DfDp_l_support
729 = periodOutArgs.supports(MEB::OUT_ARG_DfDp,period_l);
730 if (!period_DfDp_l_support.none()) {
731 outArgs.setSupports( MEB::OUT_ARG_DfDp, l, period_DfDp_l_support );
732 outArgs.set_DfDp_properties(
733 l, periodOutArgs.get_DfDp_properties(period_l) );
734 }
735 }
736
737 // DgDx_dot
738 const MEB::DerivativeSupport
739 period_DgDx_dot_support = periodOutArgs.supports(MEB::OUT_ARG_DgDx_dot,g_index_);
740 if (!period_DgDx_dot_support.none()) {
741 outArgs.setSupports( MEB::OUT_ARG_DgDx_dot, 0, period_DgDx_dot_support );
742 outArgs.set_DgDx_dot_properties(
743 0, periodOutArgs.get_DgDx_dot_properties(g_index_) );
744 }
745
746 // DgDx
747 const MEB::DerivativeSupport
748 period_DgDx_support = periodOutArgs.supports(MEB::OUT_ARG_DgDx,g_index_);
749 if (!period_DgDx_support.none()) {
750 outArgs.setSupports( MEB::OUT_ARG_DgDx, 0, period_DgDx_support );
751 outArgs.set_DgDx_properties(
752 0, periodOutArgs.get_DgDx_properties(g_index_) );
753 }
754
755 // DgDp(l)
756 for ( int l = 0; l < Np_; ++l ) {
757 const int period_l = this->period_l(l);
758 const MEB::DerivativeSupport period_DgDp_l_support
759 = periodOutArgs.supports(MEB::OUT_ARG_DgDp, g_index_, period_l);
760 if (!period_DgDp_l_support.none()) {
761 outArgs.setSupports( MEB::OUT_ARG_DgDp, 0, l, period_DgDp_l_support );
762 outArgs.set_DgDp_properties(
763 0, l, periodOutArgs.get_DgDp_properties(g_index_,period_l) );
764 }
765 }
766
767 return outArgs;
768
769}
770
771
772template<class Scalar>
773void DefaultMultiPeriodModelEvaluator<Scalar>::evalModelImpl(
774 const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
775 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
776 ) const
777{
778
779 using Teuchos::rcp_dynamic_cast;
781 typedef ModelEvaluatorBase MEB;
783
784 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
785 "DefaultMultiPeriodModelEvaluator",inArgs,outArgs,periodModel_ );
786 // ToDo: You will have to set the verbosity level for each of the
787 // periodModels_[i] individually below!
788
789 const int N = x_bar_space_->numBlocks();
790 const int Np = this->Np_;
791 //const int Ng = this->Ng_;
792
793 //
794 // A) Setup InArgs
795 //
796
797 RCP<const ProductVectorBase<Scalar> > x_bar;
798 if (inArgs.supports(MEB::IN_ARG_x)) {
799 x_bar = rcp_dynamic_cast<const ProductVectorBase<Scalar> >(
800 inArgs.get_x(), true );
802 is_null(x_bar), std::logic_error,
803 "Error, if x is supported, it must be set!"
804 );
805 }
806
807 //
808 // B) Setup OutArgs
809 //
810
811 RCP<ProductVectorBase<Scalar> > f_bar;
812 if (outArgs.supports(MEB::OUT_ARG_f)) {
813 f_bar = rcp_dynamic_cast<ProductVectorBase<Scalar> >(
814 outArgs.get_f(), true );
815 }
816
817 Array<MEB::Derivative<Scalar> > DfDp_bar(Np);
818 Array<RCP<ProductMultiVectorBase<Scalar> > > DfDp_bar_mv(Np);
819 for ( int l = 0; l < Np; ++l ) {
820 if (!outArgs.supports(MEB::OUT_ARG_DfDp,l).none()) {
821 MEB::Derivative<Scalar>
822 DfDp_bar_l = outArgs.get_DfDp(l);
823 DfDp_bar[l] = DfDp_bar_l;
824 DfDp_bar_mv[l] = rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(
825 DfDp_bar_l.getMultiVector(), true );
827 (
828 !DfDp_bar_l.isEmpty()
829 &&
830 (
831 is_null(DfDp_bar_mv[l])
832 ||
833 DfDp_bar_l.getMultiVectorOrientation() != MEB::DERIV_MV_BY_COL
834 )
835 ),
836 std::logic_error,
837 "Error, we currently can only handle DfDp as an column-based multi-vector!"
838 );
839 }
840 }
841
842 RCP<BlockedLinearOpBase<Scalar> > W_op_bar;
843 if (outArgs.supports(MEB::OUT_ARG_W_op)) {
844 W_op_bar = rcp_dynamic_cast<BlockedLinearOpBase<Scalar> >(
845 outArgs.get_W_op(), true
846 );
847 }
848
849 RCP<VectorBase<Scalar> >
850 g_bar = outArgs.get_g(0);
851
852 MEB::Derivative<Scalar> DgDx_dot_bar;
853 RCP<ProductMultiVectorBase<Scalar> > DgDx_dot_bar_mv;
854 if (!outArgs.supports(MEB::OUT_ARG_DgDx_dot,0).none()) {
855 DgDx_dot_bar = outArgs.get_DgDx_dot(0);
856 DgDx_dot_bar_mv = rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(
857 DgDx_dot_bar.getMultiVector(), true );
859 (
860 !DgDx_dot_bar.isEmpty()
861 &&
862 (
863 is_null(DgDx_dot_bar_mv)
864 ||
865 DgDx_dot_bar.getMultiVectorOrientation() != MEB::DERIV_TRANS_MV_BY_ROW
866 )
867 ),
868 std::logic_error,
869 "Error, we currently can only handle DgDx_dot as an row-based multi-vector!"
870 );
871 }
872
873 MEB::Derivative<Scalar> DgDx_bar;
874 RCP<ProductMultiVectorBase<Scalar> > DgDx_bar_mv;
875 if (!outArgs.supports(MEB::OUT_ARG_DgDx,0).none()) {
876 DgDx_bar = outArgs.get_DgDx(0);
877 DgDx_bar_mv = rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(
878 DgDx_bar.getMultiVector(), true );
880 (
881 !DgDx_bar.isEmpty()
882 &&
883 (
884 is_null(DgDx_bar_mv)
885 ||
886 DgDx_bar.getMultiVectorOrientation() != MEB::DERIV_TRANS_MV_BY_ROW
887 )
888 ),
889 std::logic_error,
890 "Error, we currently can only handle DgDx as an row-based multi-vector!"
891 );
892 }
893
894 Array<MEB::Derivative<Scalar> > DgDp_bar(Np);
895 Array<RCP<MultiVectorBase<Scalar> > > DgDp_bar_mv(Np);
896 for ( int l = 0; l < Np; ++l ) {
897 if (!outArgs.supports(MEB::OUT_ARG_DgDp,0,l).none()) {
898 MEB::Derivative<Scalar>
899 DgDp_bar_l = outArgs.get_DgDp(0,l);
900 DgDp_bar[l] = DgDp_bar_l;
901 DgDp_bar_mv[l] = DgDp_bar_l.getMultiVector();
903 !DgDp_bar_l.isEmpty() && is_null(DgDp_bar_mv[l]),
904 std::logic_error,
905 "Error, we currently can only handle DgDp as some type of multi-vector!"
906 );
907 }
908 }
909
910 //
911 // C) Evaluate the model
912 //
913
914 // C.1) Set up storage and do some initializations
915
916 MEB::InArgs<Scalar>
917 periodInArgs = periodModel_->createInArgs();
918 // ToDo: The above will have to change if you allow different structures for
919 // each period model!
920
921 // Set all of the parameters that will just be passed through
922 for ( int l = 0; l < Np; ++l )
923 periodInArgs.set_p( period_l(l), inArgs.get_p(l) ); // Can be null just fine
924
925 MEB::OutArgs<Scalar>
926 periodOutArgs = periodModel_->createOutArgs();
927 // ToDo: The above will have to change if you allow different structures for
928 // each period model!
929
930 // Create storage for period g's that will be summed into global g_bar
931 periodOutArgs.set_g(
932 g_index_, createMember<Scalar>( periodModel_->get_g_space(g_index_) ) );
933
934 // Zero out global g_bar that will be summed into below
935 if (!is_null(g_bar) )
936 assign( g_bar.ptr(), ST::zero() );
937
938 // Set up storage for peroid DgDp[l] objects that will be summed into global
939 // DgDp_bar[l] and zero out DgDp_bar[l] that will be summed into.
940 for ( int l = 0; l < Np; ++l ) {
941 if ( !is_null(DgDp_bar_mv[l]) ) {
942 assign(DgDp_bar_mv[l].ptr(), ST::zero());
943 periodOutArgs.set_DgDp(
944 g_index_, period_l(l),
945 create_DgDp_mv(
946 *periodModel_, g_index_, period_l(l),
947 DgDp_bar[l].getMultiVectorOrientation()
948 )
949 );
950 }
951 }
952
953 // C.2) Loop over periods and assemble the model
954
955 for ( int i = 0; i < N; ++i ) {
956
957 VOTSME thyraModel_outputTempState(periodModels_[i],out,verbLevel);
958
959 // C.2.a) Set period-speicific InArgs and OutArgs
960
961 for ( int k = 0; k < numPeriodZs(); ++k )
962 periodInArgs.set_p( z_indexes_[k], z_[i][k] );
963
964 if (!is_null(x_bar))
965 periodInArgs.set_x(x_bar->getVectorBlock(i));
966
967 if (!is_null(f_bar))
968 periodOutArgs.set_f(f_bar->getNonconstVectorBlock(i)); // Updated in place!
969
970 if ( !is_null(W_op_bar) )
971 periodOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,i));
972
973 for ( int l = 0; l < Np; ++l ) {
974 if ( !is_null(DfDp_bar_mv[l]) ) {
975 periodOutArgs.set_DfDp(
976 period_l(l),
977 MEB::Derivative<Scalar>(
978 DfDp_bar_mv[l]->getNonconstMultiVectorBlock(i),
979 MEB::DERIV_MV_BY_COL
980 )
981 );
982 }
983 }
984
985 if ( !is_null(DgDx_dot_bar_mv) ) {
986 periodOutArgs.set_DgDx_dot(
987 g_index_,
988 MEB::Derivative<Scalar>(
989 DgDx_dot_bar_mv->getNonconstMultiVectorBlock(i),
990 MEB::DERIV_TRANS_MV_BY_ROW
991 )
992 );
993 }
994
995 if ( !is_null(DgDx_bar_mv) ) {
996 periodOutArgs.set_DgDx(
997 g_index_,
998 MEB::Derivative<Scalar>(
999 DgDx_bar_mv->getNonconstMultiVectorBlock(i),
1000 MEB::DERIV_TRANS_MV_BY_ROW
1001 )
1002 );
1003 }
1004
1005 // C.2.b) Evaluate the period model
1006
1007 periodModels_[i]->evalModel( periodInArgs, periodOutArgs );
1008
1009 // C.2.c) Process output arguments that need processed
1010
1011 // Sum into global g_bar
1012 if (!is_null(g_bar)) {
1013 Vp_StV( g_bar.ptr(), g_weights_[i], *periodOutArgs.get_g(g_index_) );
1014 }
1015
1016 // Sum into global DgDp_bar
1017 for ( int l = 0; l < Np; ++l ) {
1018 if ( !is_null(DgDp_bar_mv[l]) ) {
1019 update(
1020 g_weights_[i],
1021 *periodOutArgs.get_DgDp(g_index_,period_l(l)).getMultiVector(),
1022 DgDp_bar_mv[l].ptr()
1023 );
1024 }
1025 }
1026
1027 // Scale DgDx_dot_bar_mv[i]
1028 if ( !is_null(DgDx_dot_bar_mv) ) {
1029 scale( g_weights_[i],
1030 DgDx_dot_bar_mv->getNonconstMultiVectorBlock(i).ptr() );
1031 }
1032
1033 // Scale DgDx_bar_mv[i]
1034 if ( !is_null(DgDx_bar_mv) ) {
1035 scale( g_weights_[i],
1036 DgDx_bar_mv->getNonconstMultiVectorBlock(i).ptr() );
1037 }
1038
1039 }
1040
1041 // ToDo: We need to do some type of global sum of g_bar and DgDp_bar to
1042 // account for other clusters of processes. I might do this with a separate
1043 // non-ANA class.
1044
1045 // Once we get here, all of the quantities should be updated and we should
1046 // be all done!
1047
1048 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
1049
1050}
1051
1052
1053// private
1054
1055
1056template<class Scalar>
1057void DefaultMultiPeriodModelEvaluator<Scalar>::set_z_indexes_and_create_period_l_map(
1058 const Array<int> &z_indexes
1059 )
1060{
1061#ifdef TEUCHOS_DEBUG
1062 TEUCHOS_TEST_FOR_EXCEPT( Np_ <= 0 && "Error, Np must be set!" );
1063#endif
1064 z_indexes_ = z_indexes;
1065 period_l_map_.resize(0);
1066 const int numTotalParams = Np_ + z_indexes_.size();
1067 Array<int>::const_iterator
1068 z_indexes_itr = z_indexes_.begin(),
1069 z_indexes_end = z_indexes_.end();
1070 int last_z_index = -1;
1071 for ( int k = 0; k < numTotalParams; ++k ) {
1072 if ( z_indexes_itr == z_indexes_end || k < *z_indexes_itr ) {
1073 // This is a "free" parameter subvector
1074 period_l_map_.push_back(k);
1075 }
1076 else {
1077 // This is a "fixed" period parameter subvector so increment
1078 // z_indexes iterator.
1079#ifdef TEUCHOS_DEBUG
1080 TEUCHOS_TEST_FOR_EXCEPT( k != *z_indexes_itr && "This should never happen!" );
1081#endif
1082 const int tmp_last_z_index = *z_indexes_itr;
1083 ++z_indexes_itr;
1084 if ( z_indexes_itr != z_indexes_end ) {
1085#ifdef TEUCHOS_DEBUG
1086 if ( last_z_index >= 0 ) {
1088 *z_indexes_itr <= last_z_index, std::logic_error,
1089 "Error, the z_indexes array = " << toString(z_indexes_)
1090 << " is not sorted or contains duplicate entries!"
1091 );
1092 }
1093#endif
1094 last_z_index = tmp_last_z_index;
1095 }
1096 }
1097 }
1098}
1099
1100
1101template<class Scalar>
1102void DefaultMultiPeriodModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
1103{
1104
1105 using Teuchos::rcp_dynamic_cast;
1106 typedef ModelEvaluatorBase MEB;
1107
1108 nominalValues_ = this->createInArgs();
1109 lowerBounds_ = this->createInArgs();
1110 upperBounds_ = this->createInArgs();
1111
1112 const MEB::InArgs<Scalar>
1113 periodNominalValues = periodModel_->getNominalValues(),
1114 periodLowerBounds = periodModel_->getLowerBounds(),
1115 periodUpperBounds = periodModel_->getUpperBounds();
1116
1117 if (periodNominalValues.supports(MEB::IN_ARG_x)) {
1118
1119 if( !is_null(periodNominalValues.get_x()) ) {
1120 // If the first peroid model has nominal values for x, then all of them
1121 // must also!
1122 RCP<Thyra::ProductVectorBase<Scalar> >
1123 x_bar_init = createProductVector(x_bar_space_);
1124 const int N = this->N();
1125 for ( int i = 0; i < N; ++i ) {
1126 assign(
1127 x_bar_init->getNonconstVectorBlock(i).ptr(),
1128 *periodModels_[i]->getNominalValues().get_x()
1129 );
1130 }
1131 nominalValues_.set_x(x_bar_init);
1132 }
1133
1134 if( !is_null(periodLowerBounds.get_x()) ) {
1135 // If the first peroid model has lower bounds for for x, then all of
1136 // them must also!
1137 RCP<Thyra::ProductVectorBase<Scalar> >
1138 x_bar_l = createProductVector(x_bar_space_);
1139 const int N = this->N();
1140 for ( int i = 0; i < N; ++i ) {
1141 assign(
1142 x_bar_l->getNonconstVectorBlock(i).ptr(),
1143 *periodModels_[i]->getLowerBounds().get_x()
1144 );
1145 }
1146 lowerBounds_.set_x(x_bar_l);
1147 }
1148
1149 if( !is_null(periodUpperBounds.get_x()) ) {
1150 // If the first peroid model has upper bounds for for x, then all of
1151 // them must also!
1152 RCP<Thyra::ProductVectorBase<Scalar> >
1153 x_bar_u = createProductVector(x_bar_space_);
1154 const int N = this->N();
1155 for ( int i = 0; i < N; ++i ) {
1156 assign(
1157 x_bar_u->getNonconstVectorBlock(i).ptr(),
1158 *periodModels_[i]->getUpperBounds().get_x()
1159 );
1160 }
1161 upperBounds_.set_x(x_bar_u);
1162 }
1163
1164 }
1165
1166 // There can only be one set of nominal values for the non-period parameters
1167 // so just take them from the first period!
1168 for ( int l = 0; l < Np_; ++l ) {
1169 const int period_l = this->period_l(l);
1170 nominalValues_.set_p(l,periodNominalValues.get_p(period_l));
1171 lowerBounds_.set_p(l,periodLowerBounds.get_p(period_l));
1172 upperBounds_.set_p(l,periodUpperBounds.get_p(period_l));
1173 }
1174
1175}
1176
1177
1178template<class Scalar>
1179RCP<ProductVectorBase<Scalar> >
1180DefaultMultiPeriodModelEvaluator<Scalar>::createProductVector(
1181 const RCP<const ProductVectorSpaceBase<Scalar> > &prodVecSpc
1182 )
1183{
1184 return Teuchos::rcp_dynamic_cast<ProductVectorBase<Scalar> >(
1185 Thyra::createMember<Scalar>(prodVecSpc), true );
1186}
1187
1188
1189} // namespace Thyra
1190
1191
1192#endif // THYRA_DEFAULT_MULTI_PERIOD_MODEL_EVALUATOR_HPP
size_type size() const
Composite subclass that takes a single ModelEvaluator object and represents it as a single aggregate ...
ArrayView< const std::string > get_g_names(int j) const
RCP< const VectorSpaceBase< Scalar > > get_x_space() const
ModelEvaluatorBase::InArgs< Scalar > getUpperBounds() const
RCP< const VectorSpaceBase< Scalar > > get_g_space(int j) const
RCP< PreconditionerBase< Scalar > > create_W_prec() const
ModelEvaluatorBase::InArgs< Scalar > getLowerBounds() const
RCP< const VectorSpaceBase< Scalar > > get_p_space(int l) const
RCP< const LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
RCP< const Array< std::string > > get_p_names(int l) const
void reset_z(const Array< Array< RCP< const VectorBase< Scalar > > > > &z)
Reset z.
RCP< const VectorSpaceBase< Scalar > > get_f_space() const
void reportFinalPoint(const ModelEvaluatorBase::InArgs< Scalar > &finalPoint, const bool wasSolved)
Ignores the final point.
void initialize(const int N, const Array< RCP< ModelEvaluator< Scalar > > > &periodModels, const Array< int > &z_indexes, const Array< Array< RCP< const VectorBase< Scalar > > > > &z, const int g_index, const Array< Scalar > g_weights, const RCP< const ProductVectorSpaceBase< Scalar > > &x_bar_space=Teuchos::null, const RCP< const ProductVectorSpaceBase< Scalar > > &f_bar_space=Teuchos::null, const RCP< LinearOpWithSolveFactoryBase< Scalar > > &W_bar_factory=Teuchos::null)
Initialize.
ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Factory interface for creating LinearOpWithSolveBase objects from compatible LinearOpBase objects.
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object.
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object.
Base subclass for ModelEvaluator that defines some basic types.
Default base class for concrete model evaluators.
Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
Abstract interface for finite-dimensional dense vectors.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool is_null(const std::shared_ptr< T > &p)
const char * toString(EConj conj)
Return a string name for a EOpTransp value. `*.
TypeTo implicit_cast(const TypeFrom &t)