00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "LOCA_MultiContinuation_ArcLengthGroup.H"
00043 #include "LOCA_MultiContinuation_ArcLengthConstraint.H"
00044 #include "LOCA_MultiContinuation_AbstractGroup.H"
00045 #include "LOCA_MultiContinuation_ConstrainedGroup.H"
00046 #include "LOCA_MultiPredictor_AbstractStrategy.H"
00047 #include "Teuchos_ParameterList.hpp"
00048 #include "LOCA_GlobalData.H"
00049 #include "NOX_Utils.H"
00050
00051 LOCA::MultiContinuation::ArcLengthGroup::ArcLengthGroup(
00052 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00053 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00054 const Teuchos::RCP<Teuchos::ParameterList>& continuationParams,
00055 const Teuchos::RCP<LOCA::MultiContinuation::AbstractGroup>& grp,
00056 const Teuchos::RCP<LOCA::MultiPredictor::AbstractStrategy>& pred,
00057 const vector<int>& paramIDs)
00058 : LOCA::MultiContinuation::ExtendedGroup(global_data, topParams,
00059 continuationParams,
00060 grp, pred, paramIDs),
00061 theta(paramIDs.size(), 1.0),
00062 doArcLengthScaling(true),
00063 gGoal(0.5),
00064 gMax(0.8),
00065 thetaMin(1.0e-3),
00066 isFirstRescale(true)
00067 {
00068 Teuchos::RCP<LOCA::MultiContinuation::ConstraintInterface> cons
00069 = Teuchos::rcp(new LOCA::MultiContinuation::ArcLengthConstraint(
00070 globalData, Teuchos::rcp(this, false)));
00071 LOCA::MultiContinuation::ExtendedGroup::setConstraints(cons, false);
00072
00073 double theta0 =
00074 continuationParams->get("Initial Scale Factor", 1.0);
00075 doArcLengthScaling =
00076 continuationParams->get("Enable Arc Length Scaling",true);
00077 gGoal =
00078 continuationParams->get("Goal Arc Length Parameter Contribution",
00079 0.5);
00080 gMax =
00081 continuationParams->get("Max Arc Length Parameter Contribution",
00082 0.8);
00083 thetaMin = continuationParams->get("Min Scale Factor", 1.0e-3);
00084
00085 for (int i=0; i<numParams; i++)
00086 theta[i] = theta0;
00087 }
00088
00089 LOCA::MultiContinuation::ArcLengthGroup::ArcLengthGroup(
00090 const LOCA::MultiContinuation::ArcLengthGroup& source,
00091 NOX::CopyType type)
00092 : LOCA::MultiContinuation::ExtendedGroup(source, type),
00093 theta(source.theta),
00094 doArcLengthScaling(source.doArcLengthScaling),
00095 gGoal(source.gGoal),
00096 gMax(source.gMax),
00097 thetaMin(source.thetaMin),
00098 isFirstRescale(source.isFirstRescale)
00099 {
00100 Teuchos::rcp_dynamic_cast<LOCA::MultiContinuation::ArcLengthConstraint>(conGroup->getConstraints())->setArcLengthGroup(Teuchos::rcp(this, false));
00101 }
00102
00103
00104 LOCA::MultiContinuation::ArcLengthGroup::~ArcLengthGroup()
00105 {
00106 }
00107
00108 NOX::Abstract::Group&
00109 LOCA::MultiContinuation::ArcLengthGroup::operator=(
00110 const NOX::Abstract::Group& source)
00111 {
00112 copy(source);
00113 return *this;
00114 }
00115
00116 Teuchos::RCP<NOX::Abstract::Group>
00117 LOCA::MultiContinuation::ArcLengthGroup::clone(NOX::CopyType type) const
00118 {
00119 return
00120 Teuchos::rcp(new LOCA::MultiContinuation::ArcLengthGroup(*this, type));
00121 }
00122
00123 void
00124 LOCA::MultiContinuation::ArcLengthGroup::copy(const NOX::Abstract::Group& src)
00125 {
00126
00127 const LOCA::MultiContinuation::ArcLengthGroup& source =
00128 dynamic_cast<const LOCA::MultiContinuation::ArcLengthGroup&>(src);
00129
00130
00131 if (this != &source) {
00132 LOCA::MultiContinuation::ExtendedGroup::copy(source);
00133 theta = source.theta;
00134 doArcLengthScaling = source.doArcLengthScaling;
00135 gGoal = source.gGoal;
00136 gMax = source.gMax;
00137 thetaMin = source.thetaMin;
00138 isFirstRescale = source.isFirstRescale;
00139 }
00140 }
00141
00142 void
00143 LOCA::MultiContinuation::ArcLengthGroup::scaleTangent()
00144 {
00145 double dpdsOld, dpdsNew;
00146 double thetaOld, thetaNew;
00147 LOCA::MultiContinuation::ExtendedVector *v, *sv;
00148
00149 scaledTangentMultiVec = tangentMultiVec;
00150
00151
00152 if (predictor->isTangentScalable()) {
00153
00154 for (int i=0; i<numParams; i++) {
00155 v =
00156 dynamic_cast<LOCA::MultiContinuation::ExtendedVector*>(&tangentMultiVec[i]);
00157 sv =
00158 dynamic_cast<LOCA::MultiContinuation::ExtendedVector*>(&scaledTangentMultiVec[i]);
00159 grpPtr->scaleVector(*(sv->getXVec()));
00160 grpPtr->scaleVector(*(sv->getXVec()));
00161
00162 if (doArcLengthScaling) {
00163
00164
00165 thetaOld = theta[i];
00166 sv->getScalars()->scale(thetaOld*thetaOld);
00167 dpdsOld = 1.0/sqrt(sv->innerProduct(*v));
00168
00169 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00170 globalData->locaUtils->out() <<
00171 std::endl << "\t" <<
00172 globalData->locaUtils->fill(64, '+') << std::endl << "\t" <<
00173 "Arc-length scaling calculation for parameter " <<
00174 getContinuationParameterName(i) << ": " << std::endl << "\t" <<
00175 "Parameter component of predictor before rescaling = " <<
00176 globalData->locaUtils->sciformat(dpdsOld) << std::endl << "\t" <<
00177 "Scale factor from previous step " <<
00178 globalData->locaUtils->sciformat(thetaOld) << std::endl << "\t" <<
00179 "Parameter contribution to arc-length equation = " <<
00180 globalData->locaUtils->sciformat(thetaOld*dpdsOld) << std::endl;
00181 }
00182
00183
00184 recalculateScaleFactor(dpdsOld, thetaOld, thetaNew);
00185
00186 sv->getScalars()->scale(thetaNew*thetaNew / (thetaOld*thetaOld));
00187
00188
00189 dpdsNew = 1.0/sqrt(sv->innerProduct(*v));
00190
00191 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00192 globalData->locaUtils->out() << std::endl << "\t" <<
00193 "Parameter component of predictor after rescaling = " <<
00194 globalData->locaUtils->sciformat(dpdsNew) << std::endl << "\t" <<
00195 "New scale factor (theta) = " <<
00196 globalData->locaUtils->sciformat(thetaNew) << std::endl << "\t" <<
00197 "Parameter contribution to arc-length equation = " <<
00198 globalData->locaUtils->sciformat(thetaNew*dpdsNew) << std::endl <<
00199 "\t" << globalData->locaUtils->fill(64, '+') << std::endl;
00200 }
00201
00202
00203 v->scale(dpdsNew);
00204 sv->scale(dpdsNew);
00205
00206 theta[i] = thetaNew;
00207
00208
00209
00210
00211
00212
00213 if (isFirstRescale) {
00214 stepSizeScaleFactor[i] = 1.0/dpdsNew;
00215 }
00216 else
00217 stepSizeScaleFactor[i] = dpdsOld/dpdsNew;
00218 }
00219 }
00220
00221 if (doArcLengthScaling && isFirstRescale)
00222 isFirstRescale = false;
00223
00224 }
00225 }
00226
00227 double
00228 LOCA::MultiContinuation::ArcLengthGroup::computeScaledDotProduct(
00229 const NOX::Abstract::Vector& x,
00230 const NOX::Abstract::Vector& y) const
00231 {
00232 const LOCA::MultiContinuation::ExtendedVector& mx =
00233 dynamic_cast<const LOCA::MultiContinuation::ExtendedVector&>(x);
00234 const LOCA::MultiContinuation::ExtendedVector& my =
00235 dynamic_cast<const LOCA::MultiContinuation::ExtendedVector&>(y);
00236
00237 double val = grpPtr->computeScaledDotProduct(*mx.getXVec(), *my.getXVec());
00238 for (int i=0; i<numParams; i++)
00239 val += theta[i] * theta[i] * mx.getScalar(i) * my.getScalar(i);
00240
00241 return val;
00242 }
00243
00244 void
00245 LOCA::MultiContinuation::ArcLengthGroup::recalculateScaleFactor(
00246 double dpds,
00247 double thetaOld,
00248 double& thetaNew)
00249 {
00250 double g = dpds*thetaOld;
00251
00252 if (g > gMax) {
00253 thetaNew = gGoal/dpds * sqrt( fabs(1.0 - g*g) / fabs(1.0 - gGoal*gGoal) );
00254
00255 if (thetaNew < thetaMin)
00256 thetaNew = thetaMin;
00257 }
00258 else
00259 thetaNew = thetaOld;
00260 }
00261
00262