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 #ifndef OPFACTORY_BINARY_ADVDIFF_SUR_HPP 45 #define OPFACTORY_BINARY_ADVDIFF_SUR_HPP 46 47 #include "ROL_Bounds.hpp" 48 #include "ROL_ScaledStdVector.hpp" 49 #include "ROL_Reduced_Objective_SimOpt.hpp" 50 #include "ROL_OptimizationProblemFactory.hpp" 51 52 #include "../../TOOLS/pdeobjective.hpp" 53 #include "../../TOOLS/pdevector.hpp" 54 #include "objsum.hpp" 55 56 template<class Real> 57 class BinaryAdvDiffFactory : public ROL::OptimizationProblemFactory<Real> { 58 private: 59 int dim_; 60 61 mutable ROL::ParameterList pl_; 62 ROL::Ptr<const Teuchos::Comm<int>> comm_; 63 ROL::Ptr<std::ostream> os_; 64 65 ROL::Ptr<FEMdata<Real>> fem_; 66 ROL::Ptr<Assembler<Real>> assembler_; 67 ROL::Ptr<ROL::Objective<Real>> penalty_, binary_; 68 ROL::Ptr<ROL::Vector<Real>> z_; 69 ROL::Ptr<ROL::BoundConstraint<Real>> bnd_; 70 71 public: BinaryAdvDiffFactory(ROL::ParameterList & pl,const ROL::Ptr<const Teuchos::Comm<int>> & comm,const ROL::Ptr<std::ostream> & os)72 BinaryAdvDiffFactory(ROL::ParameterList &pl, 73 const ROL::Ptr<const Teuchos::Comm<int>> &comm, 74 const ROL::Ptr<std::ostream> &os) 75 : pl_(pl), comm_(comm), os_(os) { 76 fem_ = ROL::makePtr<FEMdata<Real>>(comm,pl_,*os_); 77 assembler_ = fem_->getAssembler(); 78 79 std::string costType = pl_.sublist("Problem").get("Control Cost Type", "TV"); 80 std::vector<ROL::Ptr<QoI<Real>>> qoi_pen(2,ROL::nullPtr); 81 if (costType=="TV") { 82 qoi_pen[0] = ROL::makePtr<QoI_TVControl_Cost_adv_diff<Real>>(fem_->getFE2(),pl_); 83 } 84 else if (costType=="L1") { 85 qoi_pen[0] = ROL::makePtr<QoI_Control_Cost_adv_diff<Real>>(fem_->getFE2(),pl_); 86 } 87 else { 88 qoi_pen[0] = ROL::makePtr<QoI_Control_Cost_L2_adv_diff<Real>>(fem_->getFE2(),pl_); 89 } 90 qoi_pen[1] = ROL::makePtr<QoI_IntegralityControl_Cost_adv_diff<Real>>(fem_->getFE2(),pl_); 91 penalty_ = ROL::makePtr<IntegralOptObjective<Real>>(qoi_pen[0],assembler_); 92 binary_ = ROL::makePtr<IntegralOptObjective<Real>>(qoi_pen[1],assembler_); 93 // Create template control vector 94 bool usePC = pl_.sublist("Problem").get("Piecewise Constant Controls", true); 95 if (usePC) { 96 int nx = pl_.sublist("Problem").get("Number of X-Cells", 4); 97 int ny = pl_.sublist("Problem").get("Number of Y-Cells", 2); 98 Real width = pl_.sublist("Geometry").get("Width", 2.0); 99 Real height = pl_.sublist("Geometry").get("Height", 1.0); 100 Real vol = width/static_cast<Real>(nx) * height/static_cast<Real>(ny); 101 ROL::Ptr<std::vector<Real>> xvec, svec; 102 xvec = ROL::makePtr<std::vector<Real>>(nx*ny,0); 103 svec = ROL::makePtr<std::vector<Real>>(nx*ny,vol); 104 std::ifstream file; 105 std::stringstream name; 106 name << "control_" << nx << "_" << ny << ".txt"; 107 file.open(name.str()); 108 if (file.is_open()) { 109 for (int i = 0; i < nx; ++i) { 110 for (int j = 0; j < ny; ++j) { 111 file >> (*xvec)[i+j*nx]; 112 } 113 } 114 } 115 file.close(); 116 z_ = ROL::makePtr<ROL::PrimalScaledStdVector<Real>>(xvec,svec); 117 } 118 else { 119 z_ = fem_->createControlVector(pl_); 120 } 121 // Create bound constraint 122 ROL::Ptr<ROL::Vector<Real>> zlop = z_->clone(), zhip = z_->clone(); 123 zlop->setScalar(static_cast<Real>(0)); 124 zhip->setScalar(static_cast<Real>(1)); 125 bnd_ = ROL::makePtr<ROL::Bounds<Real>>(zlop,zhip); 126 } 127 update(void)128 void update(void) {} 129 buildObjective(void)130 ROL::Ptr<ROL::Objective<Real>> buildObjective(void) { 131 return ROL::makePtr<Sum_Objective<Real>>(fem_,penalty_,binary_,pl_,os_,true); 132 } 133 buildSolutionVector(void)134 ROL::Ptr<ROL::Vector<Real>> buildSolutionVector(void) { 135 ROL::Ptr<ROL::Vector<Real>> z = z_->clone(); z->set(*z_); 136 return z; 137 } 138 buildBoundConstraint(void)139 ROL::Ptr<ROL::BoundConstraint<Real>> buildBoundConstraint(void) { 140 return bnd_; 141 //ROL::Ptr<ROL::Vector<Real>> zlop = z_->clone(), zhip = z_->clone(); 142 //zlop->setScalar(static_cast<Real>(0)); 143 //zhip->setScalar(static_cast<Real>(1)); 144 //return ROL::makePtr<ROL::Bounds<Real>>(zlop,zhip); 145 //return ROL::nullPtr; 146 } 147 buildEqualityConstraint(void)148 ROL::Ptr<ROL::Constraint<Real>> buildEqualityConstraint(void) { 149 return ROL::nullPtr; 150 } 151 buildEqualityMultiplier(void)152 ROL::Ptr<ROL::Vector<Real>> buildEqualityMultiplier(void) { 153 return ROL::nullPtr; 154 } 155 buildInequalityConstraint(void)156 ROL::Ptr<ROL::Constraint<Real>> buildInequalityConstraint(void) { 157 return ROL::nullPtr; 158 } 159 buildInequalityMultiplier(void)160 ROL::Ptr<ROL::Vector<Real>> buildInequalityMultiplier(void) { 161 return ROL::nullPtr; 162 } 163 buildInequalityBoundConstraint(void)164 ROL::Ptr<ROL::BoundConstraint<Real>> buildInequalityBoundConstraint(void) { 165 return ROL::nullPtr; 166 } 167 getState(ROL::Ptr<ROL::Vector<Real>> & u,const ROL::Ptr<ROL::Vector<Real>> & z) const168 void getState(ROL::Ptr<ROL::Vector<Real>> &u, const ROL::Ptr<ROL::Vector<Real>> &z) const { 169 u = fem_->createStateVector(pl_); 170 Misfit_Objective<Real> obj(fem_,pl_); 171 obj.solvePDE(*u,*z); 172 } 173 getAssembler(void) const174 ROL::Ptr<Assembler<Real>> getAssembler(void) const { 175 return assembler_; 176 } 177 print(std::ostream & stream=std::cout)178 void print(std::ostream &stream = std::cout) { 179 assembler_->printMeshData(stream); 180 } 181 check(std::ostream & stream=std::cout)182 void check(std::ostream &stream = std::cout) { 183 ROL::Ptr<ROL::Vector<Real>> z1 = z_->clone(); z1->randomize(); 184 ROL::Ptr<ROL::Vector<Real>> z2 = z_->clone(); z2->randomize(); 185 penalty_->checkGradient(*z1,*z2,true,stream); 186 penalty_->checkHessVec(*z1,*z2,true,stream); 187 binary_->checkGradient(*z1,*z2,true,stream); 188 binary_->checkHessVec(*z1,*z2,true,stream); 189 Misfit_Objective<Real> obj(fem_,pl_); 190 obj.checkGradient(*z1,*z2,true,stream); 191 obj.checkHessVec(*z1,*z2,true,stream); 192 } 193 }; 194 195 #endif 196