continuation_system.c

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

// $Id: continuation_system.C 2956 2008-08-02 21:50:10Z benkirk $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

// LibMesh includes
#include "continuation_system.h"
#include "linear_solver.h"
#include "time_solver.h"
#include "newton_solver.h"
#include "sparse_matrix.h"

ContinuationSystem::ContinuationSystem (EquationSystems& es,
					const std::string& name,
					const unsigned int number)
  : Parent(es, name, number),
    continuation_parameter(NULL),
    quiet(true),
    continuation_parameter_tolerance(1.e-6),
    solution_tolerance(1.e-6),
    initial_newton_tolerance(0.01),
    old_continuation_parameter(0.),
    min_continuation_parameter(0.),
    max_continuation_parameter(0.),
    Theta(1.),
    Theta_LOCA(1.),
    //tau(1.),
    n_backtrack_steps(5),
    n_arclength_reductions(5),
    ds_min(1.e-8),
    predictor(Euler),
    newton_stepgrowth_aggressiveness(1.),
    newton_progress_check(true),
    rhs_mode(Residual),
    linear_solver(LinearSolver<Number>::build()),
    tangent_initialized(false),
    newton_solver(NULL),
    dlambda_ds(0.707),
    ds(0.1),
    ds_current(0.1),
    previous_dlambda_ds(0.),
    previous_ds(0.),
    newton_step(0)
{
  // Warn about using untested code
  untested();
}




ContinuationSystem::~ContinuationSystem ()
{
  this->clear();
}




void ContinuationSystem::clear()
{
  // FIXME: Do anything here, e.g. zero vectors, etc?

  // Call the Parent's clear function
  Parent::clear();
}



void ContinuationSystem::init_data ()
{
  // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
  du_ds = &(add_vector("du_ds"));

  // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
  previous_du_ds = &(add_vector("previous_du_ds"));

  // Add a vector to keep track of the previous nonlinear solution
  // at the old value of lambda.
  previous_u = &(add_vector("previous_u"));

  // Add a vector to keep track of the temporary solution "y" of Ay=G_{\lambda}.
  y = &(add_vector("y"));

  // Add a vector to keep track of the "old value" of "y" which is the solution of Ay=G_{\lambda}.
  y_old = &(add_vector("y_old"));

  // Add a vector to keep track of the temporary solution "z" of Az=-G.
  z = &(add_vector("z"));

  // Add a vector to keep track of the Newton update during the constrained PDE solves.
  delta_u = &(add_vector("delta_u"));

  // Call the Parent's initialization routine.
  Parent::init_data();
}




void ContinuationSystem::solve()
{
  // Set the Residual RHS mode, and call the normal solve routine.
  rhs_mode      = Residual;
  DifferentiableSystem::solve();
}




void ContinuationSystem::initialize_tangent()
{
  // Be sure the tangent was not already initialized.
  libmesh_assert (!tangent_initialized);
  
  // Compute delta_s_zero, the initial arclength travelled during the
  // first step.  Here we assume that previous_u and lambda_old store
  // the previous solution and control parameter.  You may need to
  // read in an old solution (or solve the non-continuation system)
  // first and call save_current_solution() before getting here.
  
  // 1.) Compute delta_s_zero as ||u|| - ||u_old|| + ...
  // Compute norms of the current and previous solutions
//   Real norm_u          = solution->l2_norm();
//   Real norm_previous_u = previous_u->l2_norm();
  
//   if (!quiet)
//     {
//       std::cout << "norm_u=" << norm_u << std::endl;
//       std::cout << "norm_previous_u=" << norm_previous_u << std::endl;
//     }
  
//   if (norm_u == norm_previous_u)
//     {
//       std::cerr << "Warning, it appears u and previous_u are the "
//   		<< "same, are you sure this is correct?"
//   		<< "It's possible you forgot to set one or the other..."
//   		<< std::endl;
//     }
  
//   Real delta_s_zero = std::sqrt(
//   				(norm_u - norm_previous_u)*(norm_u - norm_previous_u) +
//   				(*continuation_parameter-old_continuation_parameter)*
//   				(*continuation_parameter-old_continuation_parameter)
//   				);

//   // 2.) Compute delta_s_zero as ||u -u_old|| + ...
//   *delta_u = *solution;
//   delta_u->add(-1., *previous_u);
//   delta_u->close();
//   Real norm_delta_u = delta_u->l2_norm();
//   Real norm_u          = solution->l2_norm();
//   Real norm_previous_u = previous_u->l2_norm();

//   // Scale norm_delta_u by the bigger of either norm_u or norm_previous_u
//   norm_delta_u /= std::max(norm_u, norm_previous_u);
  
//   if (!quiet)
//     {
//       std::cout << "norm_u=" << norm_u << std::endl;
//       std::cout << "norm_previous_u=" << norm_previous_u << std::endl;
//       //std::cout << "norm_delta_u=" << norm_delta_u << std::endl;
//       std::cout << "norm_delta_u/max(|u|,|u_old|)=" << norm_delta_u << std::endl;
//       std::cout << "|norm_u-norm_previous_u|=" << std::abs(norm_u - norm_previous_u) << std::endl;
//     }
  
//   const Real dlambda = *continuation_parameter-old_continuation_parameter;

//   if (!quiet)
//     std::cout << "dlambda=" << dlambda << std::endl;
  
//   Real delta_s_zero = std::sqrt(
//   				(norm_delta_u*norm_delta_u) +
//   				(dlambda*dlambda)
//   				);
  
//   if (!quiet)
//     std::cout << "delta_s_zero=" << delta_s_zero << std::endl;

  // 1.) + 2.)
//   // Now approximate the initial tangent d(lambda)/ds
//   this->dlambda_ds = (*continuation_parameter-old_continuation_parameter) / delta_s_zero;


//   // We can also approximate the deriv. wrt s by finite differences:
//   // du/ds = (u1 - u0) / delta_s_zero.
//   // FIXME: Use delta_u from above if we decide to keep that method.
//   *du_ds = *solution;
//   du_ds->add(-1., *previous_u);
//   du_ds->scale(1./delta_s_zero);
//   du_ds->close();


  // 3.) Treating (u-previous_u)/(lambda - lambda_old) as an approximation to du/d(lambda),
  // we follow the same technique as Carnes and Shadid.
//   const Real dlambda = *continuation_parameter-old_continuation_parameter;
//   libmesh_assert (dlambda > 0.);
  
//   // Use delta_u for temporary calculation of du/d(lambda)
//   *delta_u = *solution;
//   delta_u->add(-1., *previous_u);
//   delta_u->scale(1. / dlambda);
//   delta_u->close();

//   // Determine initial normalization parameter
//   const Real solution_size = std::max(solution->l2_norm(), previous_u->l2_norm());
//   if (solution_size > 1.)
//     {
//       Theta = 1./solution_size;
      
//       if (!quiet)
// 	std::cout << "Setting Normalization Parameter Theta=" << Theta << std::endl;
//     }
  
//   // Compute d(lambda)/ds
//   // The correct sign of d(lambda)/ds should be positive, since we assume that (lambda > lambda_old)
//   // but we could always double-check that as well.
//   Real norm_delta_u = delta_u->l2_norm();
//   this->dlambda_ds = 1. / std::sqrt(1. + Theta*Theta*norm_delta_u*norm_delta_u);

//   // Finally, compute du/ds = d(lambda)/ds * du/d(lambda)
//   *du_ds = *delta_u;
//   du_ds->scale(dlambda_ds);
//   du_ds->close();


  // 4.) Use normalized arclength formula to estimate delta_s_zero
//   // Determine initial normalization parameter
//   set_Theta();

//   // Compute (normalized) delta_s_zero
//   *delta_u = *solution;
//   delta_u->add(-1., *previous_u);
//   delta_u->close();
//   Real norm_delta_u = delta_u->l2_norm();
  
//   const Real dlambda = *continuation_parameter-old_continuation_parameter;

//   if (!quiet)
//     std::cout << "dlambda=" << dlambda << std::endl;
  
//   Real delta_s_zero = std::sqrt(
//   				(Theta_LOCA*Theta_LOCA*Theta*norm_delta_u*norm_delta_u) +
//   				(dlambda*dlambda)
//   				);
//   *du_ds = *delta_u;
//   du_ds->scale(1./delta_s_zero);
//   dlambda_ds = dlambda / delta_s_zero;
  
//   if (!quiet)
//     {
//       std::cout << "delta_s_zero=" << delta_s_zero << std::endl;
//       std::cout << "initial d(lambda)/ds|_0 = " << dlambda_ds << std::endl;
//       std::cout << "initial ||du_ds||_0 = " << du_ds->l2_norm() << std::endl;
//     }

//   // FIXME: Also store the initial finite-differenced approximation to -du/dlambda as y.
//   // We stick to the convention of storing negative y, since that is what we typically
//   // solve for anyway.
//   *y = *delta_u;
//   y->scale(-1./dlambda);
//   y->close();

  

  // 5.) Assume dlambda/ds_0 ~ 1/sqrt(2) and determine the value of Theta_LOCA which
  // will satisfy this criterion

  // Initial change in parameter
  const Real dlambda = *continuation_parameter-old_continuation_parameter;
  libmesh_assert (dlambda != 0.0);
  
  // Ideal initial value of dlambda_ds
  dlambda_ds = 1. / std::sqrt(2.);
  if (dlambda < 0.)
    dlambda_ds *= -1.;
  
  // This also implies the initial value of ds
  ds_current = dlambda / dlambda_ds;

  if (!quiet)
    std::cout << "Setting ds_current|_0=" << ds_current << std::endl;
  
  // Set y = -du/dlambda using finite difference approximation
  *y = *solution;
  y->add(-1., *previous_u);
  y->scale(-1./dlambda);
  y->close();
  const Real ynorm=y->l2_norm();
  
  // Finally, set the value of du_ds to be used in the upcoming
  // tangent calculation. du/ds = du/dlambda * dlambda/ds
  *du_ds = *y;
  du_ds->scale(-dlambda_ds);
  du_ds->close();

  // Determine additional solution normalization parameter
  // (Since we just set du/ds, it will be:  ||du||*||du/ds||)
  set_Theta();

  // The value of Theta_LOCA which makes dlambda_ds = 1/sqrt(2),
  // assuming our Theta = ||du||^2.
  // Theta_LOCA = std::abs(dlambda);

  // Assuming general Theta
  Theta_LOCA = std::sqrt(1./Theta/ynorm/ynorm);
  
  
  if (!quiet)
    {
      std::cout << "Setting initial Theta_LOCA = " << Theta_LOCA << std::endl;
      std::cout << "Theta_LOCA^2*Theta         = " << Theta_LOCA*Theta_LOCA*Theta << std::endl;
      std::cout << "initial d(lambda)/ds|_0    = " << dlambda_ds << std::endl;
      std::cout << "initial ||du_ds||_0        = " << du_ds->l2_norm() << std::endl;
    }


  
  // OK, we estimated the tangent at point u0.
  // Now, to estimate the tangent at point u1, we call the solve_tangent routine.

  // Set the flag which tells us the method has been initialized.
  tangent_initialized = true;

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






// This is most of the "guts" of this class.  This is where we implement
// our custom Newton iterations and perform most of the solves.
void ContinuationSystem::continuation_solve()
{
  // Be sure the user has set the continuation parameter pointer
  if (!continuation_parameter)
    {
      std::cerr << "You must set the continuation_parameter pointer "