continuation_system.c

一个用来实现偏微分方程中网格的计算库

  y->close();
}



void ContinuationSystem::set_Theta()
{
  // // Use the norm of the latest solution, squared.
  //const Real normu = solution->l2_norm();
  //libmesh_assert (normu != 0.0);
  //Theta = 1./normu/normu;
  
  // // 1.) Use the norm of du, squared
//   *delta_u = *solution;
//   delta_u->add(-1, *previous_u);
//   delta_u->close();
//   const Real normdu = delta_u->l2_norm();
  
//   if (normdu < 1.) // don't divide by zero or make a huge scaling parameter.
//     Theta = 1.;
//   else
//     Theta = 1./normdu/normdu;

  // 2.) Use 1.0, i.e. don't scale
  Theta=1.;

  // 3.) Use a formula which attempts to make the "solution triangle" isosceles.
//   libmesh_assert (std::abs(dlambda_ds) < 1.);

//   *delta_u = *solution;
//   delta_u->add(-1, *previous_u);
//   delta_u->close();
//   const Real normdu = delta_u->l2_norm();

//   Theta = std::sqrt(1. - dlambda_ds*dlambda_ds) / normdu * tau * ds;


//   // 4.) Use the norm of du and the norm of du/ds
//   *delta_u = *solution;
//   delta_u->add(-1, *previous_u);
//   delta_u->close();
//   const Real normdu   = delta_u->l2_norm();
//   du_ds->close();
//   const Real normduds = du_ds->l2_norm();

//   if (normduds < 1.e-12)
//     {
//       std::cout << "Setting initial Theta= 1./normdu/normdu" << std::endl;
//       std::cout << "normdu=" << normdu << std::endl;

//       // Don't use this scaling if the solution delta is already O(1)
//       if (normdu > 1.)
// 	Theta = 1./normdu/normdu;
//       else
// 	Theta = 1.;
//     }
//   else
//     {
//       std::cout << "Setting Theta= 1./normdu/normduds" << std::endl;
//       std::cout << "normdu=" << normdu << std::endl;
//       std::cout << "normduds=" << normduds << std::endl;

//       // Don't use this scaling if the solution delta is already O(1)
//       if ((normdu>1.) || (normduds>1.))
// 	Theta = 1./normdu/normduds;
//       else
// 	Theta = 1.;
//     }
  
 if (!quiet)
   std::cout << "Setting Normalization Parameter Theta=" << Theta << std::endl;
}



void ContinuationSystem::set_Theta_LOCA()
{
  // We also recompute the LOCA normalization parameter based on the
  // most recently computed value of dlambda_ds
  // if (!quiet)
  //   std::cout << "(Theta_LOCA) dlambda_ds=" << dlambda_ds << std::endl;

  // Formula makes no sense if |dlambda_ds| > 1
  libmesh_assert (std::abs(dlambda_ds) < 1.);
  
  // 1.) Attempt to implement the method in LOCA paper
//   const Real g = 1./std::sqrt(2.); // "desired" dlambda_ds

//   // According to the LOCA people, we only renormalize for
//   // when |dlambda_ds| exceeds some pre-selected maximum (which they take to be zero, btw).
//   if (std::abs(dlambda_ds) > .9)
//     {
//       // Note the *= ... This is updating the previous value of Theta_LOCA
//       // Note: The LOCA people actually use Theta_LOCA^2 to normalize their arclength constraint.
//       Theta_LOCA *= std::abs( (dlambda_ds/g)*std::sqrt( (1.-g*g) / (1.-dlambda_ds*dlambda_ds) ) );
  
//       // Suggested max-allowable value for Theta_LOCA
//       if (Theta_LOCA > 1.e8)
// 	{
// 	  Theta_LOCA = 1.e8;

// 	  if (!quiet)
// 	    std::cout << "max Theta_LOCA=" << Theta_LOCA << " has been selected." << std::endl;
// 	}
//     }
//   else
//     Theta_LOCA=1.0;

  // 2.) FIXME: Should we do *= or just =?  This function is of dlambda_ds is
  //  < 1,  |dlambda_ds| < 1/sqrt(2) ~~ .7071
  //  > 1,  |dlambda_ds| > 1/sqrt(2) ~~ .7071
  Theta_LOCA *= std::abs( dlambda_ds / std::sqrt( (1.-dlambda_ds*dlambda_ds) ) );

  // Suggested max-allowable value for Theta_LOCA.  I've never come close
  // to this value in my code.
  if (Theta_LOCA > 1.e8)
    {
      Theta_LOCA = 1.e8;
      
      if (!quiet)
	std::cout << "max Theta_LOCA=" << Theta_LOCA << " has been selected." << std::endl;
    }
    
  // 3.) Use 1.0, i.e. don't scale
  //Theta_LOCA=1.0;

  if (!quiet)
    std::cout << "Setting Theta_LOCA=" << Theta_LOCA << std::endl;
}



void ContinuationSystem::update_solution()
{
  // Set some stream formatting flags
  unsigned int old_precision = std::cout.precision();
  std::cout.precision(16);
  std::cout.setf(std::ios_base::scientific);
  
  // We must have a tangent that makes sense before we can update the solution.
  libmesh_assert (tangent_initialized);

  // Compute tau, the stepsize scaling parameter which attempts to
  // reduce ds when the angle between the most recent two tangent
  // vectors becomes large.  tau is actually the (absolute value of
  // the) cosine of the angle between these two vectors... so if tau ~
  // 0 the angle is ~ 90 degrees, while if tau ~ 1 the angle is ~ 0
  // degrees.
  y_old->close();
  y->close();
  const Real yoldnorm = y_old->l2_norm();
  const Real ynorm = y->l2_norm();
  const Number yoldy = y_old->dot(*y);
  const Real yold_over_y = yoldnorm/ynorm;
  
  if (!quiet)
    {
      std::cout << "yoldnorm=" << yoldnorm << std::endl;
      std::cout << "ynorm="    << ynorm << std::endl;
      std::cout << "yoldy="    << yoldy << std::endl;
      std::cout << "yoldnorm/ynorm=" << yoldnorm/ynorm << std::endl;
    }      

  // Save the current value of ds before updating it
  previous_ds = ds_current;
  
  // // 1.) Cosine method (for some reason this always predicts the angle is ~0)
  // // Don't try divinding by zero
  // if ((yoldnorm > 1.e-12) && (ynorm > 1.e-12))
  //   tau = std::abs(yoldy) / yoldnorm  / ynorm;
  // else
  //   tau = 1.;

  // // 2.) Relative size of old and new du/dlambda method with cutoff of 0.9
  // if ((yold_over_y < 0.9) && (yold_over_y > 1.e-6))
  //   tau = yold_over_y;
  // else
  //   tau = 1.;

  // 3.) Grow (or shrink) the arclength stepsize by the ratio of du/dlambda, but do not
  // exceed the user-specified value of ds.
  if (yold_over_y > 1.e-6)
    {
      // // 1.) Scale current ds by the ratio of successive tangents.
      //       ds_current *= yold_over_y;
      //       if (ds_current > ds)
      // 	ds_current = ds;

      // 2.) Technique 1 tends to shrink the step fairly well (and even if it doesn't
      // get very small, we still have step reduction) but it seems to grow the step
      // very slowly.  Another possible technique is step-doubling:
//       if (yold_over_y > 1.)
//       	ds_current *= 2.;
//       else
// 	ds_current *= yold_over_y;

      // 3.) Technique 2 may over-zealous when we are also using the Newton stepgrowth
      // factor.  For technique 3 we multiply by yold_over_y unless yold_over_y > 2
      // in which case we use 2.
      //       if (yold_over_y > 2.)
      // 	ds_current *= 2.;
      //       else
      // 	ds_current *= yold_over_y;

      // 4.) Double-or-halve.  We double the arc-step if the ratio of successive tangents
      // is larger than 'double_threshold', halve it if it is less than 'halve_threshold'
      const Real double_threshold = 0.5;
      const Real halve_threshold  = 0.5;
      if (yold_over_y > double_threshold)
      	ds_current *= 2.;
      else if (yold_over_y < halve_threshold)
      	ds_current *= 0.5;
	
      
      // Also possibly use the number of Newton iterations required to compute the previous
      // step (relative to the maximum-allowed number of Newton iterations) to grow the step.
      if (newton_stepgrowth_aggressiveness > 0.)
	{
	  libmesh_assert (newton_solver != NULL);
	  const unsigned int Nmax = newton_solver->max_nonlinear_iterations;

	  // // The LOCA Newton step growth technique (note: only grows step length)
	  // const Real stepratio = static_cast<Real>(Nmax-(newton_step+1))/static_cast<Real>(Nmax-1.);
	  // const Real newtonstep_growthfactor = 1. + newton_stepgrowth_aggressiveness*stepratio*stepratio;

	  // The "Nopt/N" method, may grow or shrink the step.  Assume Nopt=Nmax/2.
	  const Real newtonstep_growthfactor =
	    newton_stepgrowth_aggressiveness * 0.5 *
	    static_cast<Real>(Nmax) / static_cast<Real>(newton_step+1);

	  if (!quiet)
	    std::cout << "newtonstep_growthfactor=" << newtonstep_growthfactor << std::endl;
	
	  ds_current *= newtonstep_growthfactor;
	}
    }

  
  // Don't let the stepsize get above the user's maximum-allowed stepsize.
  if (ds_current > ds)
    ds_current = ds;
  
  // Check also for a minimum allowed stepsize.
  if (ds_current < ds_min)
    {
      std::cout << "Enforcing minimum-allowed arclength stepsize of " << ds_min << std::endl;
      ds_current = ds_min;
    }
  
  if (!quiet)
    {
      std::cout << "Current step size: ds_current=" << ds_current << std::endl;
    }
  
  // Recompute scaling factor Theta for
  // the current solution before updating.
  set_Theta();

  // Also, recompute the LOCA scaling factor, which attempts to
  // maintain a reasonable value of dlambda/ds
  set_Theta_LOCA();
    
  std::cout << "Theta*Theta_LOCA^2=" << Theta*Theta_LOCA*Theta_LOCA << std::endl;

  // Based on the asymptotic singular behavior of du/dlambda near simple turning points,
  // we can compute a single parameter which may suggest that we are close to a singularity.
  *delta_u = *solution;
  delta_u->add(-1, *previous_u);
  delta_u->close();
  const Real normdu   = delta_u->l2_norm();
  const Real C = (std::log (Theta_LOCA*normdu) /
		  std::log (std::abs(*continuation_parameter-old_continuation_parameter))) - 1.0;
  if (!quiet)
    std::cout << "C=" << C << std::endl;
  
  // Save the current value of u and lambda before updating.
  save_current_solution();

  if (!quiet)
    {
      std::cout << "Updating the solution with the tangent guess." << std::endl;
      std::cout << "||u_old||=" << this->solution->l2_norm() << std::endl;
      std::cout << "lambda_old=" << *continuation_parameter << std::endl;
    }

  // Since we solved for the tangent vector, now we can compute an
  // initial guess for the new solution, and an initial guess for the
  // new value of lambda.
  apply_predictor();

  if (!quiet)
    {
      std::cout << "||u_new||=" << this->solution->l2_norm() << std::endl;
      std::cout << "lambda_new=" << *continuation_parameter << std::endl;
    }
  
  // Unset previous stream flags
  std::cout.precision(old_precision);
  std::cout.unsetf(std::ios_base::scientific);
}




void ContinuationSystem::save_current_solution()
{
  // Save the old solution vector
  *previous_u = *solution;

  // Save the old value of lambda
  old_continuation_parameter = *continuation_parameter;
}



void ContinuationSystem::apply_predictor()
{
  if (predictor == Euler)
    {
      // 1.) Euler Predictor
      // Predict next the solution
      solution->add(ds_current, *du_ds);
      solution->close();
  
      // Predict next parameter value
      *continuation_parameter += ds_current*dlambda_ds;
    }


  else if (predictor == AB2)
    {
      // 2.) 2nd-order explicit AB predictor
      libmesh_assert(previous_ds != 0.0);
      const Real stepratio = ds_current/previous_ds;
      
      // Build up next solution value.
      solution->add( 0.5*ds_current*(2.+stepratio), *du_ds);
      solution->add(-0.5*ds_current*stepratio     , *previous_du_ds);
      solution->close();
      
      // Next parameter value
      *continuation_parameter +=
	0.5*ds_current*((2.+stepratio)*dlambda_ds -
			stepratio*previous_dlambda_ds);
    }

  else
    {
      // Unknown predictor
      libmesh_error();
    }
  
}