Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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();
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.