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_MSG1(verb, tab0 << "=================================================="
00068 << endl
00069 << tab0 << " starting time integration " << endl
00070 << tab0 << "==================================================");
00071 PLAYA_MSG2(verb, tab1 << "Initial time: " << t);
00072 PLAYA_MSG2(verb, tab1 << "Final time: " << tStop);
00073 PLAYA_MSG2(verb, tab1 << "Initial stepsize: " << dt);
00074 PLAYA_MSG2(verb, tab1 << "Min allowed stepsize: " << minStepsize);
00075 PLAYA_MSG2(verb, tab1 << "Max allowed stepsize: " << maxStepsize);
00076 PLAYA_MSG2(verb, tab1 << "Step tolerance: " << tau);
00077 PLAYA_MSG2(verb, tab1 << "Method order: " << p);
00078 PLAYA_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_MSG2(verb, tab1 << "Terminate on event detection: "
00105 << terminateOnDetection);
00106
00107 while (t < tStop && step < maxSteps)
00108 {
00109 Tabs tab2;
00110 PLAYA_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_MSG3(verb, tab3 << "timestep exceeds interval: reducing to dt=" << dt);
00118 PLAYA_MSG3(verb, tab3 << " step " << step << " is now time=" << t
00119 << " to " << t + dt << " stepsize=" << dt);
00120 }
00121
00122 PLAYA_MSG5(verb, tab3 << "doing full step");
00123 bool stepOK = prob_.step(t, uCur, t+dt, uSolnOneFullStep, solver_);
00124 PLAYA_MSG5(verb, tab3 << "doing first half step");
00125 stepOK = prob_.step(t, uCur, t+0.5*dt, uSolnHalfStep, solver_);
00126 PLAYA_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_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_MSG1(verb, tab2 << "WARNING: minimum stepsize reached");
00145 }
00146 else
00147 {
00148 dt = safety*dtNew;
00149 }
00150 atMaxStepsize = false;
00151 PLAYA_MSG3(verb, tab3 << "reducing step to dt=" << dt);
00152 numShrink++;
00153 continue;
00154 }
00155
00156 if (eventHandler_.get())
00157 {
00158 gotEvent = eventHandler_->checkForEvent(t, uCur, t+0.5*dt, uSolnHalfStep);
00159 if (gotEvent && terminateOnDetection) break;
00160
00161 gotEvent = eventHandler_->checkForEvent(t+0.5*dt, uSolnHalfStep,
00162 t+dt, uSolnTwoHalfSteps);
00163 }
00164
00165 while ( tNextWrite >= t && tNextWrite < t+dt)
00166 {
00167 double tNode0 = t;
00168 double tNode1 = t+dt/2.0;
00169 double tNode2 = t+dt;
00170
00171 double phi0 = (tNextWrite-tNode1)*(tNextWrite-tNode2)/(tNode0-tNode1)/(tNode0-tNode2);
00172 double phi1 = (tNextWrite-tNode0)*(tNextWrite-tNode2)/(tNode1-tNode0)/(tNode1-tNode2);
00173 double phi2 = (tNextWrite-tNode0)*(tNextWrite-tNode1)/(tNode2-tNode0)/(tNode2-tNode1);
00174
00175 Vector<double> y0 = getDiscreteFunctionVector(uCur);
00176 Vector<double> y1 = getDiscreteFunctionVector(uSolnHalfStep);
00177 Vector<double> y2 = getDiscreteFunctionVector(uSolnTwoHalfSteps);
00178
00179 Vector<double> yTmp = phi0*y0 + phi1*y1 + phi2*y2;
00180 setDiscreteFunctionVector(uTmp, yTmp);
00181
00182 write(writeIndex, tNextWrite, uTmp);
00183
00184 tNextWrite += writeInterval;
00185 writeIndex++;
00186 }
00187
00188 minStepUsed = min(minStepUsed, fabs(dt));
00189 maxStepUsed = max(maxStepUsed, fabs(dt));
00190
00191 updateDiscreteFunction(uSolnTwoHalfSteps, uCur);
00192 t += dt;
00193 step++;
00194
00195 if (gotEvent && terminateOnDetection) break;
00196
00197 if (t < tStop && dtNew > dt && !atMaxStepsize)
00198 {
00199 if (dtNew >= maxStepsize)
00200 {
00201 dt = maxStepsize;
00202 atMaxStepsize = true;
00203 PLAYA_MSG2(verb, tab2 << "WARNING: maximum stepsize reached");
00204 }
00205 else
00206 {
00207 dt = dtNew;
00208 }
00209 atMinStepsize = false;
00210 PLAYA_MSG3(verb, tab3 << "increasing step to dt=" << dt);
00211 numGrow++;
00212 }
00213 }
00214
00215
00216
00217 PLAYA_MSG1(verb, tab0 << "=================================================================="
00218 << endl
00219 << tab0 << " done time integration ");
00220 PLAYA_MSG2(verb, tab0 << "------------------------------------------------------------------");
00221 PLAYA_MSG2(verb, tab1 << "Final time: " << t);
00222 PLAYA_MSG2(verb, tab1 << "Num steps: " << step);
00223 PLAYA_MSG2(verb, tab1 << "Min stepsize used: " << minStepUsed
00224 << " (would need " << fabs(t-stepControl_.tStart_)/minStepUsed
00225 << " steps of this size)");
00226 PLAYA_MSG2(verb, tab1 << "Max stepsize used: " << maxStepUsed);
00227 PLAYA_MSG2(verb, tab1 << "Num nonlinear solves: " << numSolves);
00228 PLAYA_MSG2(verb, tab1 << "Num step reductions: " << numShrink);
00229 PLAYA_MSG2(verb, tab1 << "Num step increases: " << numGrow);
00230 if (terminateOnDetection && gotEvent)
00231 {
00232 PLAYA_MSG2(verb, tab1 << "Terminated by event detection at time="
00233 << eventHandler_->eventTime());
00234 }
00235 PLAYA_MSG1(verb, tab0 << "==================================================================");
00236
00237 }
00238
00239
00240 void DoublingStepController::write(int index, double t, const Expr& u) const
00241 {
00242 Tabs tab(0);
00243 string name = outputControl_.filename_ + "-" + Teuchos::toString(index);
00244
00245 FieldWriter w = outputControl_.wf_.createWriter(name);
00246
00247 PLAYA_MSG1(outputControl_.verbosity_, tab << "writing to " << name);
00248
00249 Mesh mesh = getDiscreteFunctionMesh(u);
00250 w.addMesh(mesh);
00251 for (int i=0; i<u.size(); i++)
00252 {
00253 w.addField("output[" + Teuchos::toString(i) + "]",
00254 new ExprFieldWrapper(u[i]));
00255 }
00256 w.write();
00257
00258 }
00259
00260
00261 }
