1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) 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 lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43
44 #pragma once
45 #ifndef ROL_DYNAMICTRACKINGFEMOBJECTIVE_HPP
46 #define ROL_DYNAMICTRACKINGFEMOBJECTIVE_HPP
47
48 #include "ROL_DynamicObjective.hpp"
49 #include "ROL_VectorWorkspace.hpp"
50
51
52 /** @ingroup func_group
53 \class ROL::DynamicTrackingFEMObjective
54 \brief Defines the time-local contribution to a quadratic tracking
55 objective
56
57 \f[ f_k(u,z) = \frac{1}{2} \int\limits_{t_{k-1}}^{t_k} \| u(t)-\tilde u(t)\|^2
58 + \alpha \|z(t)\^2\,\mathm{d}t \f]
59
60 Currently approximates the state disparity norm using linear finite elements
61 which couples the old and new state contributions
62 */
63
64
65 namespace ROL {
66
67 template<typename Real>
68 class DynamicTrackingFEMObjective : public DynamicObjective<Real> {
69 public:
70
71 using V = Vector<Real>;
72 using TS = TimeStamp<Real>;
73 using size_type = typename PartitionedVector<Real>::size_type;
74
75 private:
76
77 Ptr<PartitionedVector<Real>> target_;
78
79 size_type Nt_; // Number of time steps
80 Real alpha_; // Regularization parameter
81
82 mutable VectorWorkspace<Real> workspace_;
83
84 public:
85
DynamicTrackingFEMObjective(const Ptr<PartitionedVector<Real>> & target,Real alpha=0.0)86 DynamicTrackingFEMObjective( const Ptr<PartitionedVector<Real>>& target, Real alpha=0.0 ) :
87 target_(target), Nt_(target_->numVectors()), alpha_(alpha) {}
88
~DynamicTrackingFEMObjective()89 virtual ~DynamicTrackingFEMObjective() {}
90
value(const V & uo,const V & un,const V & z,const TS & timeStamp) const91 virtual Real value( const V& uo, const V& un,
92 const V& z, const TS& timeStamp ) const override {
93
94 Real dt = timeStamp.t.at(1)-timeStamp.t.at(0);
95 Real w = 2.0*dt/3.0;
96
97 size_type k = timeStamp.k;
98
99 auto res_n = workspace_.copy(un);
100 auto res_o = workspace_.copy(uo);
101
102 Real result = 0.5*dt*alpha_*z.dot(z);
103
104 res_n->set(un);
105 res_n->axpy( -1.0, *(target_->get(k)) );
106 result += w*res_n->dot(*res_n);
107
108 if( k>0 ) {
109 res_o->set(uo);
110 res_o->axpy( -1, *(target_->get(k-1) ) );
111 result += w*res_n->dot(*res_o);
112 result += w*res_o->dot(*res_o);
113 }
114
115 return result;
116 }
117
118 //----------------------------------------------------------------------------
119 // Gradient Terms
gradient_uo(V & g,const V & uo,const V & un,const V & z,const TS & timeStamp) const120 virtual void gradient_uo( V& g, const V& uo, const V& un,
121 const V& z, const TS& timeStamp ) const override {
122 Real dt = timeStamp.t.at(1)-timeStamp.t.at(0);
123 Real w = dt/3.0;
124
125 size_type k = timeStamp.k;
126
127 auto res_n = workspace_.copy(un);
128 auto res_o = workspace_.copy(uo);
129
130
131
132
133 else g.zero();
134 }
135
gradient_un(V & g,const V & uo,const V & un,const V & z,const TS & timeStamp) const136 virtual void gradient_un( V& g, const V& uo, const V& un,
137 const V& z, const TS& timeStamp ) const override {
138 Real dt = timeStamp.t.at(1)-timeStamp.t.at(0);
139 g.set(un);
140 g.axpy(-1.0, *(target_->get(timeStamp.k)) );
141 g.scale(0.5*dt);
142 }
143
gradient_z(V & g,const V & uo,const V & un,const V & z,const TS & timeStamp) const144 virtual void gradient_z( V& g, const V& uo, const V& un,
145 const V& z, const TS& timeStamp ) const override {
146 Real dt = timeStamp.t.at(1)-timeStamp.t.at(0);
147 g.set(z);
148 g.scale(dt*alpha_);
149 }
150
151 //----------------------------------------------------------------------------
152 // Hessian-Vector product terms
hessVec_uo_uo(V & hv,const V & v,const V & uo,const V & un,const V & z,const TS & timeStamp) const153 virtual void hessVec_uo_uo( V& hv, const V& v, const V& uo, const V& un,
154 const V& z, const TS& timeStamp ) const override {
155 Real dt = timeStamp.t.at(1)-timeStamp.t.at(0);
156 if( timeStamp.k>0 ) {
157 hv.set(v);
158 hv.axpy(-1.0, *(target_->get(timeStamp.k-1)) );
159 hv.scale(0.5*dt);
160 }
161 else hv.zero();
162 }
163
hessVec_uo_un(V & hv,const V & v,const V & uo,const V & un,const V & z,const TS & timeStamp) const164 virtual void hessVec_uo_un( V& hv, const V& v, const V& uo, const V& un,
165 const V& z, const TS& timeStamp ) const override {
166 hv.zero();
167 }
168
hessVec_uo_z(V & hv,const V & v,const V & uo,const V & un,const V & z,const TS & timeStamp) const169 virtual void hessVec_uo_z( V& hv, const V& v, const V& uo, const V& un,
170 const V& z, const TS& timeStamp ) const override {
171 hv.zero();
172 }
173
hessVec_un_uo(V & hv,const V & v,const V & uo,const V & un,const V & z,const TS & timeStamp) const174 virtual void hessVec_un_uo( V& hv, const V& v, const V& uo, const V& un,
175 const V& z, const TS& timeStamp ) const override {
176 hv.zero();
177 }
178
hessVec_un_un(V & hv,const V & v,const V & uo,const V & un,const V & z,const TS & timeStamp) const179 virtual void hessVec_un_un( V& hv, const V& v, const V& uo, const V& un,
180 const V& z, const TS& timeStamp ) const override {
181 hv.set(v);
182 hv.scale(0.5*(timeStamp.t.at(1)-timeStamp.t.at(0)));
183 }
184
hessVec_un_z(V & hv,const V & v,const V & uo,const V & un,const V & z,const TS & timeStamp) const185 virtual void hessVec_un_z( V& hv, const V& v, const V& uo, const V& un,
186 const V& z, const TS& timeStamp ) const override {
187 hv.zero();
188 }
189
hessVec_z_uo(V & hv,const V & v,const V & uo,const V & un,const V & z,const TS & timeStamp) const190 virtual void hessVec_z_uo( V& hv, const V& v, const V& uo, const V& un,
191 const V& z, const TS& timeStamp ) const override {
192 hv.zero();
193 }
194
hessVec_z_un(V & hv,const V & v,const V & uo,const V & un,const V & z,const TS & timeStamp) const195 virtual void hessVec_z_un( V& hv, const V& v, const V& uo, const V& un,
196 const V& z, const TS& timeStamp ) const override {
197 hv.zero();
198 }
199
hessVec_z_z(V & hv,const V & v,const V & uo,const V & un,const V & z,const TS & timeStamp) const200 virtual void hessVec_z_z( V& hv, const V& v, const V& uo, const V& un,
201 const V& z, const TS& timeStamp ) const override {
202 hv.set(v);
203 hv.scale(alpha_*(timeStamp.t.at(1)-timeStamp.t.at(0)));
204 }
205
206 }; // DynamicTrackingFEMObjective
207
208
209 template<typename Real>
210 inline Ptr<DynamicObjective<Real>>
make_DynamicTrackingFEMObjective(const Ptr<PartitionedVector<Real>> & target,Real alpha=0.0)211 make_DynamicTrackingFEMObjective( const Ptr<PartitionedVector<Real>>& target, Real alpha=0.0 ) {
212 Ptr<DynamicObjective<Real>> obj = makePtr<DynamicTrackingFEMObjective<Real>>(target,alpha);
213 return obj;
214 }
215
216 } // namespace ROL
217
218 #endif // ROL_DYNAMICTRACKINGFEMOBJECTIVE_HPP
219
220
221