1 /*! 2 * \file mfront/src/PowellDogLegAlgorithmBase.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \brief 22 août 2014 6 * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights 7 * reserved. 8 * This project is publicly released under either the GNU GPL Licence 9 * or the CECILL-A licence. A copy of thoses licences are delivered 10 * with the sources of TFEL. CEA or EDF may also distribute this 11 * project under specific licensing conditions. 12 */ 13 14 #include<iostream> 15 #include<sstream> 16 17 #include"TFEL/Raise.hxx" 18 #include"TFEL/Utilities/StringAlgorithms.hxx" 19 #include"MFront/BehaviourDescription.hxx" 20 #include"MFront/NonLinearSystemSolverBase.hxx" 21 #include"MFront/PowellDogLegAlgorithmBase.hxx" 22 23 namespace mfront{ 24 25 std::vector<std::string> getReservedNames()26 PowellDogLegAlgorithmBase::getReservedNames() 27 { 28 return {"pdl_g","pdl_g2","pdl_0", 29 "pdl_1","pdl_2","pdl_3", 30 "pdl_alpha","pdl_cste", 31 "powell_dogleg_trust_region_size_inv"}; 32 } // end of PowellDogLegAlgorithmBase::getReservedNames 33 34 std::pair<bool,PowellDogLegAlgorithmBase::tokens_iterator> treatSpecificKeywords(BehaviourDescription & mb,const std::string & key,const tokens_iterator p,const tokens_iterator pe)35 PowellDogLegAlgorithmBase::treatSpecificKeywords(BehaviourDescription& mb, 36 const std::string& key, 37 const tokens_iterator p, 38 const tokens_iterator pe) 39 { 40 using namespace tfel::material; 41 using namespace tfel::utilities; 42 const auto h = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 43 if(key=="@PowellDogLegTrustRegionSize"){ 44 CxxTokenizer::TokensContainer::const_iterator current = p; 45 CxxTokenizer::checkNotEndOfLine("PowellDogLegAlgorithmBase::treatSpecificKeywords",current,pe); 46 CxxTokenizer::checkNotEndOfLine("ImplicitDSLBase::treatPowellDogLegTrustRegionSize" 47 "Cannot read pdl_trs value.",current,pe); 48 const auto pdl_trs = convert<double>(current->value); 49 tfel::raise_if(pdl_trs<0,"ImplicitDSLBase::treatPowellDogLegTrustRegionSize: " 50 "region size must be positive."); 51 ++current; 52 CxxTokenizer::readSpecifiedToken("ImplicitDSLBase::treatPowellDogLegTrustRegionSize",";",current,pe); 53 mb.addParameter(h,VariableDescription("real","powell_dogleg_trust_region_size",1u,0u)); 54 mb.setParameterDefaultValue(h,"powell_dogleg_trust_region_size",pdl_trs); 55 return {true,current}; 56 } 57 return {false,p}; 58 } // end of PowellDogLegAlgorithmBase::treatSpecificKeywords 59 60 completeVariableDeclaration(BehaviourDescription & mb)61 void PowellDogLegAlgorithmBase::completeVariableDeclaration(BehaviourDescription& mb) 62 { 63 using namespace tfel::material; 64 const auto h = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 65 if(!mb.hasParameter(h,"powell_dogleg_trust_region_size")){ 66 mb.addParameter(h,VariableDescription("real","powell_dogleg_trust_region_size",1u,0u)); 67 mb.setParameterDefaultValue(h,"powell_dogleg_trust_region_size",1.e-4); 68 } 69 } // end of PowellDogLegAlgorithmBase::completeVariableDeclaration 70 71 void writePowellDogLegStep(std::ostream & out,const BehaviourDescription & mb,const Hypothesis h,const std::string & B,const std::string & f,const std::string & pn)72 PowellDogLegAlgorithmBase::writePowellDogLegStep(std::ostream& out, 73 const BehaviourDescription& mb, 74 const Hypothesis h, 75 const std::string& B, 76 const std::string& f, 77 const std::string& pn) 78 { 79 const auto& d = mb.getBehaviourData(h); 80 const auto n = d.getIntegrationVariables().getTypeSize(); 81 out << "if(abs(" << pn<< ")<(" << n << ")*(this->powell_dogleg_trust_region_size)){\n" 82 << "// using the newton method only\n"; 83 NonLinearSystemSolverBase::writeLimitsOnIncrementValues(out,mb,h,pn); 84 out << "this->zeros -= " << pn<< ";\n" 85 << "} else { \n" 86 << "// computing the steepest descent step\n" 87 << "tvector<" << n << ",real> pdl_g;\n" 88 << "tvector<" << n << ",real> pdl_g2;\n" 89 << "for(unsigned short idx=0;idx!=" << n << ";++idx){\n" 90 << "pdl_g[idx]=real(0);\n" 91 << "for(unsigned short idx2=0;idx2!=" << n << ";++idx2){\n" 92 << "pdl_g[idx] += (" << B << "(idx2,idx)) * (" << f << "(idx2));\n" 93 << "}\n" 94 << "}\n" 95 << "for(unsigned short idx=0;idx!=" << n << ";++idx){\n" 96 << "pdl_g2[idx]=real(0);\n" 97 << "for(unsigned short idx2=0;idx2!=" << n << ";++idx2){\n" 98 << "pdl_g2[idx] += (" << B << "(idx,idx2)) * pdl_g(idx2);\n" 99 << "}\n" 100 << "}\n" 101 << "const real pdl_cste = (pdl_g|pdl_g)/(pdl_g2|pdl_g2);\n" 102 << "pdl_g *= pdl_cste;\n" 103 << "if(abs(pdl_g)<(" << n << ")*(this->powell_dogleg_trust_region_size)){\n" 104 << "const real pdl_0 = (this->powell_dogleg_trust_region_size)*(this->powell_dogleg_trust_region_size);\n" 105 << "const real pdl_1 = (pdl_g|pdl_g);\n" 106 << "const real pdl_2 = ((" << pn << ")|pdl_g);\n" 107 << "const real pdl_3 = ((" << pn << ")|(" << pn << "));\n" 108 << "const real pdl_alpha = " 109 << "(pdl_0-pdl_1)/((pdl_2-pdl_1)+sqrt(max((pdl_2-pdl_0)*(pdl_2-pdl_0)+(pdl_3-pdl_0)*(pdl_0-pdl_1),real(0))));\n" 110 << "pdl_g = pdl_alpha*(" << pn<< ") + (1-pdl_alpha)*pdl_g;\n" 111 << "} else {\n" 112 << "const real pdl_alpha = (this->powell_dogleg_trust_region_size)/(norm(pdl_g));\n" 113 << "pdl_g *= pdl_alpha;\n" 114 << "}\n"; 115 NonLinearSystemSolverBase::writeLimitsOnIncrementValues(out,mb,h,"pdl_g"); 116 out << "this->zeros -= pdl_g;\n" 117 << "}\n"; 118 } // end of ImplicitDSLBase::writePowellDogLegStep 119 120 PowellDogLegAlgorithmBase::~PowellDogLegAlgorithmBase() = default; 121 122 } // end of namespace mfront 123