adaptive unsteady newton.cpp

double Problem::
adaptive_unsteady_newton_solve(const double &dt_desired,
                               const double &epsilon,
                               const bool &shift_values)
{
 //First, we need to backup the existing dofs, in case the timestep is
 //rejected

 //Find total number of dofs on current processor
 unsigned n_dof_local = dof_distribution_pt()->nrow_local();

 // Backup current time step dofs and time in case timestep is rejected.
 Vector<double> dofs_current(n_dof_local);
 for(unsigned i=0;i<n_dof_local;i++) dofs_current[i] = dof(i);
 double time_current = time_pt()->time();

 //Flag to detect whether the timestep has been rejected or not
 unsigned REJECT_TIMESTEP=0;
 //Flag to detect whether any of the timesteppers are adaptive
 unsigned ADAPTIVE_FLAG=0;

 //Determine the number of timesteppers
 unsigned n_time_steppers = ntime_stepper();
 //Find out whether any of the timesteppers are adaptive
 for(unsigned i=0;i<n_time_steppers;i++)
  {
   if(time_stepper_pt(i)->adaptive_flag())
    {
     ADAPTIVE_FLAG=1;
     break;
    }
  }

 // Our own copy of dt_desired that we can modify
 double dt_actual = dt_desired;

 //This loop surrounds the adaptive time-stepping critera
 do
  {
   // Initially we assume that this step will succeed
   REJECT_TIMESTEP=0;

   // and that this dt value is ok
   double dt_rescaling_factor = 1.0;

   //Attempt to solve the non-linear system
   try
    {
     //Solve the non-linear problem at this timestep
     unsteady_newton_solve(dt_actual, false, ADAPTIVE_FLAG);
    }
   //Catch any exceptions thrown
   catch(NewtonSolverError &error)
    {
     //If it's a solver error then die
     if(error.linear_solver_error)
      { //??ds should have its own error class?
       std::string error_message =
        "USER-DEFINED ERROR IN NEWTON SOLVER\n";
       error_message +=  "ERROR IN THE LINEAR SOLVER\n";

       //Die
       throw OomphLibError(error_message,
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
     else
      {
       //Reject the timestep, if we have an exception
       oomph_info << "TIMESTEP REJECTED" << std::endl;
       REJECT_TIMESTEP=1;
         
       // Half the time step
       dt_rescaling_factor = 0.5;
      }
    }

   // If we have an adapative timestepper (and we haven't already failed)
   // then calculate the error estimate and rescaling factor.
   if(ADAPTIVE_FLAG && !REJECT_TIMESTEP)
    {
     //Once timestep has been accepted can do fancy error processing
     //Set the error weights
     for(unsigned i=0;i<n_time_steppers;i++)
      {
       time_stepper_pt(i)->set_error_weights();
      }

     // Get a global error norm to use in adaptivity (as specified by the
     // problem sub-class writer). Prevent a divide by zero if the solution
     // gives very close to zero error. Error norm should never be negative
     // but use absolute value just in case.
     double error = std::max(std::abs(global_temporal_error_norm()),
                             1e-12);

     //Calculate the scaling  factor
     dt_rescaling_factor = 
      std::pow((epsilon/error), (1.0/(1.0+time_stepper_pt()->order())));

     oomph_info << "Timestep scaling factor is  " << dt_rescaling_factor << std::endl;
     oomph_info << "Estimated timestepping error is " << error << std::endl;
    } //End of if adaptive flag



   // Calculate the next time step size and check it's ok
   // ============================================================
   
   // Calculate the possible next time step, if no error conditions
   // trigger.
   double new_dt_candidate = dt_rescaling_factor * dt_actual;

   // Check that the scaling factor is within the allowed range
   if(dt_rescaling_factor > DTSF_max_increase)
    {
     oomph_info << "Tried to increase dt by the ratio " << dt_rescaling_factor
                << " which is above the maximum (" << DTSF_max_increase
                << "). Attempting to increase by the maximum ratio instead.";
     new_dt_candidate = DTSF_max_increase * dt_actual;
    }
   // If we have already rejected the timestep then don't do this check
   // because DTSF will definitely be too small.
   else if((!REJECT_TIMESTEP) && (dt_rescaling_factor <= DTSF_min_decrease))
    {
     oomph_info << "Timestep would decrease by " << dt_rescaling_factor
                << " which is less than the minimum scaling factor " 
                << DTSF_min_decrease << std::endl
                << "TIMESTEP REJECTED" << std::endl;
     REJECT_TIMESTEP=1;
    }

   // Now check that the new dt is within the allowed range
   if(new_dt_candidate > Maximum_dt)
    {
     oomph_info << "Tried to increase dt to " << new_dt_candidate
                << " which is above the maximum (" << Maximum_dt
                << "). I increased it to the maximum value instead.";
     dt_actual = Maximum_dt;
    }
   else if(new_dt_candidate < Minimum_dt)
    {
     std::ostringstream err;
     err << "Tried to reduce dt to " << new_dt_candidate 
         << " which is less than the minimum dt (" << Minimum_dt
         << ")." << std::endl;
     throw OomphLibError(err.str(), OOMPH_EXCEPTION_LOCATION,
                         OOMPH_CURRENT_FUNCTION);
    }
   else
    {
     dt_actual = new_dt_candidate;
    }


   // If we are rejecting this attempt then revert the dofs etc.
   if(REJECT_TIMESTEP)
    {
     //Reset the time
     time_pt()->time() = time_current;

     //Reload the dofs
     unsigned ni = dofs_current.size();
     for(unsigned i=0; i<ni; i++) {dof(i) = dofs_current[i];}

#ifdef OOMPH_HAS_MPI
     // Synchronise the solution on different processors (on each submesh)
     this->synchronise_all_dofs();
#endif

     //Call all "after" actions, e.g. to handle mesh updates
     actions_after_newton_step();
     actions_before_newton_convergence_check();
     actions_after_newton_solve();
     actions_after_implicit_timestep();
    }

  }
 //Keep this loop going until we accept the timestep
 while(REJECT_TIMESTEP);

 // Once the timestep has been accepted, return the time step that should be
 // used next time.
 return dt_actual;
}

void Problem::unsteady_newton_solve(const double &dt, const bool &shift_values,
                                    const bool& calculate_predictions)
{
 //Shift the time values and the dts, according to the control flag
 if(shift_values) {shift_time_values();}

 // Advance global time and set current value of dt
 time_pt()->time()+=dt;
 time_pt()->dt()=dt;

 //Find out how many timesteppers there are
 unsigned n_time_steppers = ntime_stepper();

 //Loop over them all and set the weights
 for(unsigned i=0;i<n_time_steppers;i++)
  {
   time_stepper_pt(i)->set_weights();
  }

 // If requested then set up predictor weights and calculate the predicted
 // values and positions.
 if(calculate_predictions)
  {
   for(unsigned i=0;i<n_time_steppers;i++)
    {
     time_stepper_pt(i)->set_predictor_weights();
    }
   this->calculate_predictions();
  }

 //Now update anything that needs updating before the timestep
 //This could be time-dependent boundary conditions, for example.
 actions_before_implicit_timestep();

 //Solve the non-linear problem for this timestep with Newton's method
 newton_solve();

 //Now update anything that needs updating after the timestep
 actions_after_implicit_timestep();
}