continuation_system.c

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

	  // step is where you "just were" and the current residual is where
	  // you are now.  It can occur that ||du||/||R|| < 1, but these are
	  // likely not good cases to attempt backtracking (?).
	  const Real norm_du_norm_R = delta_u->l2_norm() / nonlinear_residual_afterstep;
	  if (!quiet)
	    std::cout << "  norm_du_norm_R=" << norm_du_norm_R << std::endl;
      
      
	  // Factor to decrease the stepsize by for backtracking
	  Real newton_stepfactor = 1.;

	  const bool attempt_backtracking =
	    (nonlinear_residual_afterstep > solution_tolerance)
	    && (nonlinear_residual_afterstep > nonlinear_residual_beforestep)
	    && (n_backtrack_steps>0)
	    && (norm_du_norm_R > 1.)
	    ;
    
	  // If residual is not reduced, do Newton back tracking.  
	  if (attempt_backtracking)
	    {
	      if (!quiet)
		std::cout << "Newton step did not reduce residual." << std::endl;

	      // back off the previous step.
	      solution->add(-1., *delta_u);
	      solution->close();
	      *continuation_parameter -= delta_lambda;

	      // Backtracking: start cutting the Newton stepsize by halves until
	      // the new residual is actually smaller...
	      for (unsigned int backtrack_step=0; backtrack_step<n_backtrack_steps; ++backtrack_step)
		{
		  newton_stepfactor *= 0.5;

		  if (!quiet)
		    std::cout << "Shrinking step size by " << newton_stepfactor << std::endl;
	      
		  // Take fractional step
		  solution->add(newton_stepfactor, *delta_u);
		  solution->close();
		  *continuation_parameter += newton_stepfactor*delta_lambda;
	      
		  rhs_mode = Residual;
		  assembly(true,   // Residual
			   false); // Jacobian
		  rhs->close();
		  nonlinear_residual_afterstep = rhs->l2_norm();

		  if (!quiet)
		    std::cout << "At shrink step "
			      << backtrack_step
			      << ", nonlinear_residual_afterstep="
			      << nonlinear_residual_afterstep
			      << std::endl;

		  if (nonlinear_residual_afterstep < nonlinear_residual_beforestep)
		    {
		      if (!quiet)
			std::cout << "Backtracking succeeded!" << std::endl;
		      
		      break; // out of backtracking loop
		    }
		  
		  else
		    {
		      // Back off that step
		      solution->add(-newton_stepfactor, *delta_u);
		      solution->close();
		      *continuation_parameter -= newton_stepfactor*delta_lambda;
		    }
	      
		  // Save a copy of the solution from before the Newton step.
		  //AutoPtr<NumericVector<Number> > prior_iterate = solution->clone();
		}
	    } // end if (attempte_backtracking)


	  // If we tried backtracking but the residual is still not reduced, print message.
	  if ((attempt_backtracking) && (nonlinear_residual_afterstep > nonlinear_residual_beforestep))
	    {
	      //std::cerr << "Backtracking failed." << std::endl;
	      std::cout << "Backtracking failed." << std::endl;
	  
	      // 1.) Quit, exit program.
	      //libmesh_error();

	      // 2.) Continue with last newton_stepfactor
	      if (newton_step<3)
		{
		  solution->add(newton_stepfactor, *delta_u);
		  solution->close();
		  *continuation_parameter += newton_stepfactor*delta_lambda;
		  if (!quiet)
		    std::cout << "Backtracking could not reduce residual ... continuing anyway!" << std::endl;
		}
	      
	      // 3.) Break out of Newton iteration loop with newton_converged = false,
	      //     reduce the arclength stepsize, and try again.
	      else
		{
		  break; // out of Newton iteration loop, with newton_converged=false
		}
	    }

	  // Another type of convergence check: suppose the residual has not been reduced
	  // from its initial value after half of the allowed Newton steps have occurred.
	  // In our experience, this typically means that it isn't going to converge and
	  // we could probably save time by dropping out of the Newton iteration loop and
	  // trying a smaller arcstep.
	  if (this->newton_progress_check)
	    {
	      if ((nonlinear_residual_afterstep > nonlinear_residual_firststep) &&
		  (newton_step+1 > static_cast<unsigned int>(0.5*newton_solver->max_nonlinear_iterations)))
		{
		  std::cout << "Progress check failed: the current residual: "
			    << nonlinear_residual_afterstep
			    << ", is\n"
			    << "larger than the initial residual, and half of the allowed\n"
			    << "number of Newton iterations have elapsed.\n"
			    << "Exiting Newton iterations with converged==false." << std::endl;
	      
		  break; // out of Newton iteration loop, newton_converged = false
		}
	    }
	  
	  // Safety check: Check the current continuation parameter against user-provided min-allowable parameter value
	  if (*continuation_parameter < min_continuation_parameter)
	    {
	      std::cout << "Continuation parameter fell below min-allowable value." << std::endl;
	      // libmesh_error();
	      break; // out of Newton iteration loop, newton_converged = false
	    }
	  
	  // Safety check: Check the current continuation parameter against user-provided max-allowable parameter value
	  if ( (max_continuation_parameter != 0.0) &&
	       (*continuation_parameter > max_continuation_parameter) )
	    {
	      std::cout << "Current continuation parameter value: "
			<< *continuation_parameter
			<< " exceeded max-allowable value."
			<< std::endl;
	      // libmesh_error();
	      break; // out of Newton iteration loop, newton_converged = false
	    }
	  
      
	  // Check the convergence of the parameter and the solution.  If they are small
	  // enough, we can break out of the Newton iteration loop.  
	  const Real norm_delta_u = delta_u->l2_norm();
	  const Real norm_u = solution->l2_norm();
	  std::cout << "  delta_lambda                   = " << delta_lambda << std::endl;
	  std::cout << "  newton_stepfactor*delta_lambda = " << newton_stepfactor*delta_lambda << std::endl;
	  std::cout << "  lambda_current                 = " << *continuation_parameter << std::endl;
	  std::cout << "  ||delta_u||                    = " << norm_delta_u << std::endl;
	  std::cout << "  ||delta_u||/||u||              = " << norm_delta_u / norm_u << std::endl;
      

	  // Evaluate the residual at the current Newton iterate.  We don't want to detect
	  // convergence due to a small Newton step when the residual is still not small.
	  rhs_mode = Residual;
	  assembly(true,   // Residual
		   false); // Jacobian
	  rhs->close();
	  const Real norm_residual = rhs->l2_norm();
	  std::cout << "  ||R||_{L2} = " << norm_residual << std::endl;
	  std::cout << "  ||R||_{L2}/||u|| = " << norm_residual / norm_u << std::endl;
      
      
	  // FIXME: The norm_delta_u tolerance (at least) should be relative.
	  // It doesn't make sense to converge a solution whose size is ~ 10^5 to
	  // a tolerance of 1.e-6.  Oh, and we should also probably check the
	  // (relative) size of the residual as well, instead of just the step.
	  if ((std::abs(delta_lambda) < continuation_parameter_tolerance) &&
	      //(norm_delta_u       < solution_tolerance)               && // This is a *very* strict criterion we can probably skip
	      (norm_residual      < solution_tolerance))
	    {
	      if (!quiet)
		std::cout << "Newton iterations converged!" << std::endl;

	      newton_converged = true;
	      break; // out of Newton iterations
	    }
	} // end nonlinear loop

      if (!newton_converged)
	{
	  std::cout << "Newton iterations of augmented system did not converge!" << std::endl;
	  
	  // Reduce ds_current, recompute the solution and parameter, and continue to next
	  // arcstep, if there is one.
	  ds_current *= 0.5;

	  // Go back to previous solution and parameter value.
	  *solution = *previous_u;
	  *continuation_parameter = old_continuation_parameter;

	  // Compute new predictor with smaller ds
	  apply_predictor();
	}
      else
	{
	  // Set step convergence and break out
	  arcstep_converged=true;
	  break; // out of arclength reduction loop
	}
      
    } // end loop over arclength reductions

  // Check for convergence of the whole arcstep.  If not converged at this
  // point, we have no choice but to quit.
  if (!arcstep_converged)
    {
      std::cout << "Arcstep failed to converge after max number of reductions! Exiting..." << std::endl;
      libmesh_error();
    }
  
  // Print converged solution control parameter and max value.
  std::cout << "lambda_current=" << *continuation_parameter << std::endl;
  //std::cout << "u_max=" << solution->max() << std::endl;

  // Reset old stream precision and flags.
  std::cout.precision(old_precision);
  std::cout.unsetf(std::ios_base::scientific);

  // Note: we don't want to go on to the next guess yet, since the user may
  // want to post-process this data.  It's up to the user to call advance_arcstep()
  // when they are ready to go on.
}



void ContinuationSystem::advance_arcstep()
{
  // Solve for the updated tangent du1/ds, d(lambda1)/ds
  solve_tangent();

  // Advance the solution and the parameter to the next value.
  update_solution();
}



// This function solves the tangent system:
// [ G_u                G_{lambda}        ][(du/ds)_new      ] = [  0 ]
// [ Theta*(du/ds)_old  (dlambda/ds)_old  ][(dlambda/ds)_new ]   [-N_s]
// The solution is given by:
// .) Let G_u y = G_lambda, then 
// .) 2nd row yields:
//    (dlambda/ds)_new = 1.0 / ( (dlambda/ds)_old - Theta*(du/ds)_old*y )
// .) 1st row yields
//    (du_ds)_new = -(dlambda/ds)_new * y
void ContinuationSystem::solve_tangent()
{
  // We shouldn't call this unless the current tangent already makes sense.
  libmesh_assert (tangent_initialized);
  
  // Set pointer to underlying Newton solver
  if (!newton_solver)
    newton_solver =
      dynamic_cast<NewtonSolver*> (this->time_solver->diff_solver().get());
  libmesh_assert (newton_solver != NULL);

  // Assemble the system matrix AND rhs, with rhs = G_{\lambda}
  this->rhs_mode = G_Lambda;
  
  // Assemble Residual and Jacobian
  this->assembly(true,   // Residual
		 true); // Jacobian

  // Not sure if this is really necessary
  rhs->close();
      
  // Solve G_u*y =  G_{\lambda}
  std::pair<unsigned int, Real> rval = 
    linear_solver->solve(*matrix,
			 *y,
			 *rhs,
			 1.e-12, // relative linear tolerance
			 2*newton_solver->max_linear_iterations);   // max linear iterations

  // FIXME: If this doesn't converge at all, the new tangent vector is
  // going to be really bad...
  
  if (!quiet)
    std::cout << "G_u*y = G_{lambda} solver converged at step "
	      << rval.first
	      << " linear tolerance = "
	      << rval.second
	      << "."
	      << std::endl;

  // Save old solution and parameter tangents for possible use in higher-order
  // predictor schemes.
  previous_dlambda_ds = dlambda_ds;
  *previous_du_ds     = *du_ds;

  
  // 1.) Previous, probably wrong, technique!
//   // Solve for the updated d(lambda)/ds
//   // denom = N_{lambda}   - (du_ds)^t y
//   //       = d(lambda)/ds - (du_ds)^t y
//   Real denom = dlambda_ds - du_ds->dot(*y);

//   //std::cout << "denom=" << denom << std::endl;
//   libmesh_assert (denom != 0.0);
  
//   dlambda_ds = 1.0 / denom;


//   if (!quiet)
//     std::cout << "dlambda_ds=" << dlambda_ds << std::endl;
  
//   // Compute the updated value of du/ds = -_dlambda_ds * y
//   du_ds->zero();
//   du_ds->add(-dlambda_ds, *y);
//   du_ds->close();

  
  // 2.) From Brian Carnes' paper...
  // According to Carnes, y comes from solving G_u * y = -G_{\lambda}
  y->scale(-1.);
  const Real ynorm = y->l2_norm();
  dlambda_ds = 1. / std::sqrt(1. + Theta_LOCA*Theta_LOCA*Theta*ynorm*ynorm);

  // Determine the correct sign for dlambda_ds.
  
  // We will use delta_u to temporarily compute this sign.
  *delta_u = *solution;
  delta_u->add(-1., *previous_u);
  delta_u->close();

  const Real sgn_dlambda_ds =
    libmesh_real(Theta_LOCA*Theta_LOCA*Theta*y->dot(*delta_u) +
    (*continuation_parameter-old_continuation_parameter));

  if (sgn_dlambda_ds < 0.)
    {
      if (!quiet)
	std::cout << "dlambda_ds is negative." << std::endl;
      
      dlambda_ds *= -1.;
    }

  // Finally, set the new tangent vector, du/ds = dlambda/ds * y.
  du_ds->zero();
  du_ds->add(dlambda_ds, *y);
  du_ds->close();

  if (!quiet)
    {
      std::cout << "d(lambda)/ds = " << dlambda_ds << std::endl;
      std::cout << "||du_ds||    = " << du_ds->l2_norm() << std::endl;
    }

  // Our next solve expects y ~ -du/dlambda, so scale it back by -1 again now.
  y->scale(-1.);