continuation_system.c

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

		<< "to a member variable of the derived class, preferably in the "
		<< "Derived class's init_data function.  This is how the ContinuationSystem "
		<< "updates the continuation parameter."
		<< std::endl;
      
      libmesh_error();
    }

  // Use extra precision for all the numbers printed in this function.
  unsigned int old_precision = std::cout.precision();
  std::cout.precision(16);
  std::cout.setf(std::ios_base::scientific);
  
  // We can't start solving the augmented PDE system unless the tangent
  // vectors have been initialized.  This only needs to occur once.
  if (!tangent_initialized)
    initialize_tangent();

  // Save the old value of -du/dlambda.  This will be used after the Newton iterations
  // to compute the angle between previous tangent vectors.  This cosine of this angle is
  //
  // tau := abs( (du/d(lambda)_i , du/d(lambda)_{i-1}) / (||du/d(lambda)_i|| * ||du/d(lambda)_{i-1}||) )
  //
  // The scaling factor tau (which should vary between 0 and 1) is used to shrink the step-size ds
  // when we are approaching a turning point.  Note that it can only shrink the step size.  
  *y_old = *y;
  
  // 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);
  
  // A pair for catching return values from linear system solves.
  std::pair<unsigned int, Real> rval;

  // Convergence flag for the entire arcstep
  bool arcstep_converged = false;

  // Begin loop over arcstep reductions.
  for (unsigned int ns=0; ns<n_arclength_reductions; ++ns)
    {
      if (!quiet)
	{
	  std::cout << "Current arclength stepsize, ds_current=" << ds_current << std::endl;
	  std::cout << "Current parameter value, lambda=" << *continuation_parameter << std::endl;
	}
      
      // Upon exit from the nonlinear loop, the newton_converged flag
      // will tell us the convergence status of Newton's method.
      bool newton_converged = false;

      // The nonlinear residual before *any* nonlinear steps have been taken.
      Real nonlinear_residual_firststep = 0.;

      // The nonlinear residual from the current "k" Newton step, before the Newton step
      Real nonlinear_residual_beforestep = 0.;

      // The nonlinear residual from the current "k" Newton step, after the Newton step
      Real nonlinear_residual_afterstep = 0.;

      // The linear solver tolerance, can be updated dynamically at each Newton step.
      Real current_linear_tolerance = 0.;
      
      // The nonlinear loop
      for (newton_step=0; newton_step<newton_solver->max_nonlinear_iterations; ++newton_step)
	{
	  std::cout << "\n === Starting Newton step " << newton_step << " ===" << std::endl;
	  
	  // Set the linear system solver tolerance
// 	  // 1.) Set the current linear tolerance based as a multiple of the current residual of the system.
// 	  const Real residual_multiple = 1.e-4;
// 	  Real current_linear_tolerance = residual_multiple*nonlinear_residual_beforestep;

// 	  // But if the current residual isn't small, don't let the solver exit with zero iterations!
// 	  if (current_linear_tolerance > 1.)
// 	    current_linear_tolerance = residual_multiple;

	  // 2.) Set the current linear tolerance based on the method based on technique of Eisenstat & Walker.
	  if (newton_step==0)
	    {
	      // At first step, only try reducing the residual by a small amount
	      current_linear_tolerance = initial_newton_tolerance;//0.01;
	    }

	  else
	    {
	      // The new tolerance is based on the ratio of the most recent tolerances
	      const Real alp=0.5*(1.+std::sqrt(5.));
	      const Real gam=0.9;

	      libmesh_assert (nonlinear_residual_beforestep != 0.0);
	      libmesh_assert (nonlinear_residual_afterstep != 0.0);

	      current_linear_tolerance = std::min(gam*std::pow(nonlinear_residual_afterstep/nonlinear_residual_beforestep, alp),
						  current_linear_tolerance*current_linear_tolerance
						  );

	      // Don't let it get ridiculously small!!
	      if (current_linear_tolerance < 1.e-12)
		current_linear_tolerance = 1.e-12;
	    }

	  if (!quiet)
	    std::cout << "Using current_linear_tolerance=" << current_linear_tolerance << std::endl;

	  
	  // Assemble the residual (and Jacobian).
	  rhs_mode = Residual;
	  assembly(true,   // Residual
		   true); // Jacobian
	  rhs->close();

	  // Save the current nonlinear residual.  We don't need to recompute the residual unless
	  // this is the first step, since it was already computed as part of the convergence check
	  // at the end of the last loop iteration.
	  if (newton_step==0)
	    {
	      nonlinear_residual_beforestep = rhs->l2_norm();

	      // Store the residual before any steps have been taken.  This will *not*
	      // be updated at each step, and can be used to see if any progress has
	      // been made from the initial residual at later steps.
	      nonlinear_residual_firststep = nonlinear_residual_beforestep;

	      const Real old_norm_u = solution->l2_norm();
	      std::cout << "  (before step) ||R||_{L2} = " << nonlinear_residual_beforestep << std::endl;
	      std::cout << "  (before step) ||R||_{L2}/||u|| = " << nonlinear_residual_beforestep / old_norm_u << std::endl;

	      // In rare cases (very small arcsteps), it's possible that the residual is
	      // already below our absolute linear tolerance.
	      if (nonlinear_residual_beforestep  < solution_tolerance)
		{
		  if (!quiet)
		    std::cout << "Initial guess satisfied linear tolerance, exiting with zero Newton iterations!" << std::endl;

		  // Since we go straight from here to the solve of the next tangent, we
		  // have to close the matrix before it can be assembled again.
		  matrix->close();
		  newton_converged=true;
		  break; // out of Newton iterations, with newton_converged=true
		}
	    }

	  else
	    {
	      nonlinear_residual_beforestep = nonlinear_residual_afterstep;
	    }

	  
	  // Solve the linear system G_u*z = G
	  // Initial guess?
	  z->zero(); // It seems to be extremely important to zero z here, otherwise the solver quits early.
	  z->close();
	  
	  // It's possible that we have selected the current_linear_tolerance so large that
	  // a guess of z=zero yields a linear system residual |Az + R| small enough that the
	  // linear solver exits in zero iterations.  If this happens, we will reduce the
	  // current_linear_tolerance until the linear solver does at least 1 iteration.
	  do
	    {
	      rval =
		linear_solver->solve(*matrix,
				     *z,
				     *rhs,
				     //1.e-12,
				     current_linear_tolerance,
				     newton_solver->max_linear_iterations);   // max linear iterations
	      
	      if (rval.first==0)
		{
		  if (newton_step==0)
		    {
		      std::cout << "Repeating initial solve with smaller linear tolerance!" << std::endl;
		      current_linear_tolerance *= initial_newton_tolerance; // reduce the linear tolerance to force the solver to do some work
		    }
		  else
		    {
		      // We shouldn't get here ... it means the linear solver did no work on a Newton
		      // step other than the first one.  If this happens, we need to think more about our
		      // tolerance selection.
		      libmesh_error();
		    }
		}
	      
	    } while (rval.first==0);

	  
	  if (!quiet)
	    std::cout << "  G_u*z = G solver converged at step "
		      << rval.first
		      << " linear tolerance = "
		      << rval.second
		      << "."
		      << std::endl;

	  // Sometimes (I am not sure why) the linear solver exits after zero iterations.
	  // Perhaps it is hitting PETSc's divergence tolerance dtol???  If this occurs,
	  // we should break out of the Newton iteration loop because nothing further is
	  // going to happen...  Of course if the tolerance is already small enough after
	  // zero iterations (how can this happen?!) we should not quit.
	  if ((rval.first == 0) && (rval.second > current_linear_tolerance*nonlinear_residual_beforestep))
	    {
	      if (!quiet)
		std::cout << "Linear solver exited in zero iterations!" << std::endl;

	      // Try to find out the reason for convergence/divergence
	      linear_solver->print_converged_reason();
	      
	      break; // out of Newton iterations
	    }
	  
	  // Note: need to scale z by -1 since our code always solves Jx=R
	  // instead of Jx=-R.
	  z->scale(-1.);
	  z->close();





	  
	  // Assemble the G_Lambda vector, skip residual.
	  rhs_mode = G_Lambda;

	  // Assemble both rhs and Jacobian
	  assembly(true,  // Residual
		   false); // Jacobian

	  // Not sure if this is really necessary
	  rhs->close();
	  const Real yrhsnorm=rhs->l2_norm();
	  if (yrhsnorm == 0.0)
	    {
	      std::cout << "||G_Lambda|| = 0" << std::endl;
	      libmesh_error();
	    }

	  // We select a tolerance for the y-system which is based on the inexact Newton
	  // tolerance but scaled by an extra term proportional to the RHS (which is not -> 0 in this case)
	  const Real ysystemtol=current_linear_tolerance*(nonlinear_residual_beforestep/yrhsnorm);
	  if (!quiet)
	    std::cout << "ysystemtol=" << ysystemtol << std::endl;
	  
	  // Solve G_u*y = G_{\lambda}
	  // FIXME: Initial guess?  This is really a solve for -du/dlambda so we could try
	  // initializing it with the latest approximation to that... du/dlambda ~ du/ds * ds/dlambda
	  //*y = *solution;
	  //y->add(-1., *previous_u);
	  //y->scale(-1. / (*continuation_parameter - old_continuation_parameter)); // Be careful of divide by zero...
	  //y->close();
	  
	  //	  const unsigned int max_attempts=1;
	  // unsigned int attempt=0;
	  // 	  do
	  // 	    {
	  // 	      if (!quiet)
	  // 		std::cout << "Trying to solve tangent system, attempt " << attempt << std::endl;
	      
	      rval =
		linear_solver->solve(*matrix,
				     *y,
				     *rhs,
				     //1.e-12, 
				     ysystemtol,
				     newton_solver->max_linear_iterations);   // max linear iterations

	      if (!quiet)
		std::cout << "  G_u*y = G_{lambda} solver converged at step "
			  << rval.first
			  << ", linear tolerance = "
			  << rval.second
			  << "."
			  << std::endl;

	      // Sometimes (I am not sure why) the linear solver exits after zero iterations.
	      // Perhaps it is hitting PETSc's divergence tolerance dtol???  If this occurs,
	      // we should break out of the Newton iteration loop because nothing further is
	      // going to happen...
	      if ((rval.first == 0) && (rval.second > ysystemtol))
		{
		  if (!quiet)
		    std::cout << "Linear solver exited in zero iterations!" << std::endl;

		  break; // out of Newton iterations
		}
	      
// 	      ++attempt;
// 	    } while ((attempt<max_attempts) && (rval.first==newton_solver->max_linear_iterations));
	    
      
      
      
      
	  // Compute N, the residual of the arclength constraint eqn.
	  // Note 1: N(u,lambda,s) := (u-u_{old}, du_ds) + (lambda-lambda_{old}, dlambda_ds) - _ds 
	  // We temporarily use the delta_u vector as a temporary vector for this calculation.
	  *delta_u = *solution;
	  delta_u->add(-1., *previous_u);

	  // First part of the arclength constraint
	  const Number N1 = Theta_LOCA*Theta_LOCA*Theta*delta_u->dot(*du_ds);
	  const Number N2 = ((*continuation_parameter) - old_continuation_parameter)*dlambda_ds;
	  const Number N3 = ds_current;

	  if (!quiet)
	    {
	      std::cout << "  N1=" << N1 << std::endl;
	      std::cout << "  N2=" << N2 << std::endl;
	      std::cout << "  N3=" << N3 << std::endl;
	    }
      
	  // The arclength constraint value 
	  const Number N = N1+N2-N3;
      
	  if (!quiet)
	    std::cout << "  N=" << N << std::endl;

	  const Number duds_dot_z = du_ds->dot(*z);
	  const Number duds_dot_y = du_ds->dot(*y);

	  //std::cout << "duds_dot_z=" << duds_dot_z << std::endl;
	  //std::cout << "duds_dot_y=" << duds_dot_y << std::endl;
	  //std::cout << "dlambda_ds=" << dlambda_ds << std::endl;

	  const Number delta_lambda_numerator   = -(N          + Theta_LOCA*Theta_LOCA*Theta*duds_dot_z);
	  const Number delta_lambda_denominator =  (dlambda_ds - Theta_LOCA*Theta_LOCA*Theta*duds_dot_y);

	  libmesh_assert (delta_lambda_denominator != 0.0);

	  // Now, we are ready to compute the step delta_lambda
	  const Number delta_lambda_comp = delta_lambda_numerator /
                                           delta_lambda_denominator;
          // Lambda is real-valued
          const Real delta_lambda = libmesh_real(delta_lambda_comp);

	  // Knowing delta_lambda, we are ready to update delta_u
	  // delta_u = z - delta_lambda*y
	  delta_u->zero();
	  delta_u->add(1., *z);
	  delta_u->add(-delta_lambda, *y);
	  delta_u->close();

	  // Update the system solution and the continuation parameter.
	  solution->add(1., *delta_u);
	  solution->close();
	  *continuation_parameter += delta_lambda;

	  // Did the Newton step actually reduce the residual?
	  rhs_mode = Residual;
	  assembly(true,   // Residual
		   false); // Jacobian
	  rhs->close();
	  nonlinear_residual_afterstep = rhs->l2_norm();

	  
	  // In a "normal" Newton step, ||du||/||R|| > 1 since the most recent