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 #include "SundanceDoublingStepController.hpp"
00032 #include "SundanceDiscreteFunction.hpp"
00033 #include "SundanceExprFieldWrapper.hpp"
00034 #include "PlayaLinearCombinationImpl.hpp"
00035 #include "PlayaOut.hpp"
00036 #include "PlayaTabs.hpp"
00037
00038 namespace Sundance
00039 {
00040
00041
00042 bool DoublingStepController::run() const
00043 {
00044 Tabs tab0(0);
00045 Tabs tab1;
00046 int verb = stepControl_.verbosity_;
00047 if (MPIComm::world().getRank()!=0) verb = 0;
00048 double t = stepControl_.tStart_;
00049 double p = stepControl_.stepOrder_;
00050 double tau = stepControl_.tau_;
00051 double tStop = stepControl_.tStop_;
00052 int maxSteps = stepControl_.maxSteps_;
00053
00054 double safety = stepControl_.stepsizeReductionSafetyFactor_;
00055
00056
00057 double dt = stepControl_.initialStepsize_;
00058
00059 double minStepsize = dt * stepControl_.minStepsizeFactor_;
00060 double maxStepsize = dt * stepControl_.maxStepsizeFactor_;
00061 bool atMinStepsize = false;
00062 bool atMaxStepsize = false;
00063
00064 double minStepUsed = fabs(dt);
00065 double maxStepUsed = fabs(dt);
00066
00067 PLAYA_ROOT_MSG1(verb, tab0 << "=================================================="
00068 << endl
00069 << tab0 << " starting time integration " << endl
00070 << tab0 << "==================================================");
00071 PLAYA_ROOT_MSG2(verb, tab1 << "Initial time: " << t);
00072 PLAYA_ROOT_MSG2(verb, tab1 << "Final time: " << tStop);
00073 PLAYA_ROOT_MSG2(verb, tab1 << "Initial stepsize: " << dt);
00074 PLAYA_ROOT_MSG2(verb, tab1 << "Min allowed stepsize: " << minStepsize);
00075 PLAYA_ROOT_MSG2(verb, tab1 << "Max allowed stepsize: " << maxStepsize);
00076 PLAYA_ROOT_MSG2(verb, tab1 << "Step tolerance: " << tau);
00077 PLAYA_ROOT_MSG2(verb, tab1 << "Method order: " << p);
00078 PLAYA_ROOT_MSG2(verb, tab1 << "Max steps: " << maxSteps);
00079
00080 Expr uCur = copyDiscreteFunction(prob_.uCur());
00081 Expr uTmp = copyDiscreteFunction(uCur);
00082 Expr uSolnOneFullStep = copyDiscreteFunction(uCur);
00083 Expr uSolnHalfStep = copyDiscreteFunction(uCur);
00084 Expr uSolnTwoHalfSteps = copyDiscreteFunction(uCur);
00085
00086 int writeIndex = 1;
00087 double tNextWrite = outputControl_.writeInterval_;
00088 int writeVerb = outputControl_.verbosity_;
00089 double writeInterval = outputControl_.writeInterval_;
00090
00091
00092 write(0, 0.0, uCur);
00093
00094
00095 int step = 0;
00096 int numShrink = 0;
00097 int numGrow = 0;
00098 int numSolves = 0;
00099
00100 bool gotEvent = false;
00101 bool terminateOnDetection = false;
00102 if (eventHandler_.get()) terminateOnDetection
00103 = eventHandler_->terminateOnDetection();
00104 PLAYA_ROOT_MSG2(verb, tab1 << "Terminate on event detection: "
00105 << terminateOnDetection);
00106
00107 while (t < tStop && step < maxSteps)
00108 {
00109 Tabs tab2;
00110 PLAYA_ROOT_MSG3(verb, tab2 << " step " << step << " time=" << t
00111 << " to " << t + dt << " stepsize=" << dt);
00112 Tabs tab3;
00113
00114 if (t + dt > tStop)
00115 {
00116 dt = tStop - t;
00117 PLAYA_ROOT_MSG3(verb, tab3 << "timestep exceeds interval: reducing to dt=" << dt);
00118 PLAYA_ROOT_MSG3(verb, tab3 << " step " << step << " is now time=" << t
00119 << " to " << t + dt << " stepsize=" << dt);
00120 }
00121
00122 PLAYA_ROOT_MSG5(verb, tab3 << "doing full step");
00123 bool stepOK = prob_.step(t, uCur, t+dt, uSolnOneFullStep, solver_);
00124 PLAYA_ROOT_MSG5(verb, tab3 << "doing first half step");
00125 stepOK = prob_.step(t, uCur, t+0.5*dt, uSolnHalfStep, solver_);
00126 PLAYA_ROOT_MSG5(verb, tab3 << "doing second half step");
00127 stepOK = prob_.step(t+0.5*dt, uSolnHalfStep, t+dt, uSolnTwoHalfSteps, solver_);
00128 numSolves += 3;
00129
00130 double err = compare_->diff(uSolnTwoHalfSteps, uSolnOneFullStep);
00131
00132 PLAYA_ROOT_MSG3(verb, tab3 << "error=" << err << " tau=" << tau);
00133
00134 double dtNew = dt;
00135 if (err < tau) dtNew = dt*pow(tau/err, 1.0/(p+1.0));
00136 else dtNew = dt*pow(tau/err, 1.0/p);
00137
00138 if (dtNew < dt && !atMinStepsize)
00139 {
00140 if (safety*dtNew <= minStepsize)
00141 {
00142 dt = minStepsize;
00143 atMinStepsize = true;
00144 PLAYA_ROOT_MSG1(verb, tab2 << "WARNING: minimum stepsize reached");
00145 }
00146 else
00147 {
00148 dt = safety*dtNew;
00149 }
00150 atMaxStepsize = false;
00151 PLAYA_ROOT_MSG3(verb, tab3 << "reducing step to dt=" << dt);
00152 numShrink++;
00153 continue;
00154 }
00155
00156 if (stepHook_.get()) stepHook_->call(t+dt, uSolnTwoHalfSteps);
00157
00158 if (eventHandler_.get())
00159 {
00160 gotEvent = eventHandler_->checkForEvent(t, uCur, t+0.5*dt, uSolnHalfStep);
00161 if (gotEvent && terminateOnDetection) break;
00162
00163 gotEvent = eventHandler_->checkForEvent(t+0.5*dt, uSolnHalfStep,
00164 t+dt, uSolnTwoHalfSteps);
00165 }
00166
00167 while ( tNextWrite >= t && tNextWrite < t+dt)
00168 {
00169 double tNode0 = t;
00170 double tNode1 = t+dt/2.0;
00171 double tNode2 = t+dt;
00172
00173 double phi0 = (tNextWrite-tNode1)*(tNextWrite-tNode2)/(tNode0-tNode1)/(tNode0-tNode2);
00174 double phi1 = (tNextWrite-tNode0)*(tNextWrite-tNode2)/(tNode1-tNode0)/(tNode1-tNode2);
00175 double phi2 = (tNextWrite-tNode0)*(tNextWrite-tNode1)/(tNode2-tNode0)/(tNode2-tNode1);
00176
00177 Vector<double> y0 = getDiscreteFunctionVector(uCur);
00178 Vector<double> y1 = getDiscreteFunctionVector(uSolnHalfStep);
00179 Vector<double> y2 = getDiscreteFunctionVector(uSolnTwoHalfSteps);
00180
00181 Vector<double> yTmp = phi0*y0 + phi1*y1 + phi2*y2;
00182 setDiscreteFunctionVector(uTmp, yTmp);
00183
00184 write(writeIndex, tNextWrite, uTmp);
00185
00186 tNextWrite += writeInterval;
00187 writeIndex++;
00188 }
00189
00190 minStepUsed = min(minStepUsed, fabs(dt));
00191 maxStepUsed = max(maxStepUsed, fabs(dt));
00192
00193 updateDiscreteFunction(uSolnTwoHalfSteps, uCur);
00194 t += dt;
00195 step++;
00196
00197 if (gotEvent && terminateOnDetection) break;
00198
00199 if (t < tStop && dtNew > dt && !atMaxStepsize)
00200 {
00201 if (dtNew >= maxStepsize)
00202 {
00203 dt = maxStepsize;
00204 atMaxStepsize = true;
00205 PLAYA_ROOT_MSG2(verb, tab2 << "WARNING: maximum stepsize reached");
00206 }
00207 else
00208 {
00209 dt = dtNew;
00210 }
00211 atMinStepsize = false;
00212 PLAYA_ROOT_MSG3(verb, tab3 << "increasing step to dt=" << dt);
00213 numGrow++;
00214 }
00215 }
00216
00217
00218
00219 PLAYA_ROOT_MSG1(verb, tab0 << "=================================================================="
00220 << endl
00221 << tab0 << " done time integration ");
00222 PLAYA_ROOT_MSG2(verb, tab0 << "------------------------------------------------------------------");
00223 PLAYA_ROOT_MSG2(verb, tab1 << "Final time: " << t);
00224 PLAYA_ROOT_MSG2(verb, tab1 << "Num steps: " << step);
00225 PLAYA_ROOT_MSG2(verb, tab1 << "Min stepsize used: " << minStepUsed
00226 << " (would need " << fabs(t-stepControl_.tStart_)/minStepUsed
00227 << " steps of this size)");
00228 PLAYA_ROOT_MSG2(verb, tab1 << "Max stepsize used: " << maxStepUsed);
00229 PLAYA_ROOT_MSG2(verb, tab1 << "Num nonlinear solves: " << numSolves);
00230 PLAYA_ROOT_MSG2(verb, tab1 << "Num step reductions: " << numShrink);
00231 PLAYA_ROOT_MSG2(verb, tab1 << "Num step increases: " << numGrow);
00232 if (terminateOnDetection && gotEvent)
00233 {
00234 PLAYA_ROOT_MSG2(verb, tab1 << "Terminated by event detection at time="
00235 << eventHandler_->eventTime());
00236 }
00237 PLAYA_ROOT_MSG1(verb, tab0 << "==================================================================");
00238
00239 }
00240
00241
00242 void DoublingStepController::write(int index, double t, const Expr& u) const
00243 {
00244 Tabs tab(0);
00245 string name = outputControl_.filename_ + "-" + Teuchos::toString(index);
00246
00247 FieldWriter w = outputControl_.wf_.createWriter(name);
00248
00249 PLAYA_ROOT_MSG1(outputControl_.verbosity_, tab << "writing to " << name);
00250
00251 Mesh mesh = getDiscreteFunctionMesh(u);
00252 w.addMesh(mesh);
00253 for (int i=0; i<u.size(); i++)
00254 {
00255 w.addField("output[" + Teuchos::toString(i) + "]",
00256 new ExprFieldWrapper(u[i]));
00257 }
00258 w.write();
00259
00260 }
00261
00262
00263 }
