libMesh::EulerSolver Class Reference

#include <euler_solver.h>

Inheritance diagram for libMesh::EulerSolver:

List of all members.

Public Types

typedef UnsteadySolver Parent
typedef DifferentiableSystem sys_type

Public Member Functions

 EulerSolver (sys_type &s)
virtual ~EulerSolver ()
virtual Real error_order () const
virtual bool element_residual (bool request_jacobian, DiffContext &)
virtual bool side_residual (bool request_jacobian, DiffContext &)
virtual void init ()
virtual void init_data ()
virtual void reinit ()
virtual void solve ()
virtual void advance_timestep ()
virtual void adjoint_advance_timestep ()
virtual void retrieve_timestep ()
Number old_nonlinear_solution (const dof_id_type global_dof_number) const
virtual Real du (const SystemNorm &norm) const
virtual bool is_steady () const
virtual void before_timestep ()
const sys_typesystem () const
sys_typesystem ()
virtual AutoPtr< DiffSolver > & diff_solver ()
virtual AutoPtr< LinearSolver
< Number > > & 
linear_solver ()
void set_solution_history (const SolutionHistory &_solution_history)
bool is_adjoint () const
void set_is_adjoint (bool _is_adjoint_value)

Static Public Member Functions

static std::string get_info ()
static void print_info (std::ostream &out=libMesh::out)
static unsigned int n_objects ()
static void enable_print_counter_info ()
static void disable_print_counter_info ()

Public Attributes

Real theta
AutoPtr< NumericVector< Number > > old_local_nonlinear_solution
bool quiet
unsigned int reduce_deltat_on_diffsolver_failure

Protected Types

typedef std::map< std::string,
std::pair< unsigned int,
unsigned int > > 
Counts

Protected Member Functions

void increment_constructor_count (const std::string &name)
void increment_destructor_count (const std::string &name)

Protected Attributes

bool first_solve
bool first_adjoint_step
AutoPtr< DiffSolver_diff_solver
AutoPtr< LinearSolver< Number > > _linear_solver
sys_type_system
AutoPtr< SolutionHistorysolution_history

Static Protected Attributes

static Counts _counts
static Threads::atomic
< unsigned int > 
_n_objects
static Threads::spin_mutex _mutex
static bool _enable_print_counter = true

Detailed Description

This class defines a theta-method Euler (defaulting to Backward Euler with theta = 1.0) solver to handle time integration of DifferentiableSystems.

This class is part of the new DifferentiableSystem framework, which is still experimental. Users of this framework should beware of bugs and future API changes.

Author:
Roy H. Stogner 2006

Definition at line 45 of file euler_solver.h.


Member Typedef Documentation

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts [protected, inherited]

Data structure to log the information. The log is identified by the class name.

Definition at line 113 of file reference_counter.h.

The parent class

Definition at line 51 of file euler_solver.h.

The type of system

Reimplemented in libMesh::EigenTimeSolver, and libMesh::SteadySolver.

Definition at line 66 of file time_solver.h.


Constructor & Destructor Documentation

libMesh::EulerSolver::EulerSolver ( sys_type s  )  [explicit]

Constructor. Requires a reference to the system to be solved.

Definition at line 27 of file euler_solver.C.

00028  : UnsteadySolver(s), theta(1.)
00029 {
00030 }

libMesh::EulerSolver::~EulerSolver (  )  [virtual]

Destructor.

Definition at line 34 of file euler_solver.C.

00035 {
00036 }


Member Function Documentation

void libMesh::UnsteadySolver::adjoint_advance_timestep (  )  [virtual, inherited]

This method advances the adjoint solution to the previous timestep, after an adjoint_solve() has been performed. This will be done before every UnsteadySolver::adjoint_solve().

Reimplemented from libMesh::TimeSolver.

Definition at line 169 of file unsteady_solver.C.

References libMesh::TimeSolver::_system, libMesh::DifferentiableSystem::deltat, libMesh::UnsteadySolver::first_adjoint_step, libMesh::System::get_dof_map(), libMesh::DofMap::get_send_list(), libMesh::System::get_vector(), libMesh::NumericVector< T >::localize(), libMesh::UnsteadySolver::old_local_nonlinear_solution, libMesh::TimeSolver::solution_history, and libMesh::System::time.

00170 {
00171   // On the first call of this function, we dont save the adjoint solution or
00172   // decrement the time, we just call the retrieve function below
00173   if(!first_adjoint_step)
00174     {
00175       // Call the store function to store the last adjoint before decrementing the time
00176       solution_history->store();
00177       // Decrement the system time
00178       _system.time -= _system.deltat;
00179     }
00180   else
00181     {
00182       first_adjoint_step = false;
00183     }
00184 
00185   // Retrieve the primal solution vectors at this time using the
00186   // solution_history object
00187   solution_history->retrieve();
00188 
00189   // Dont forget to localize the old_nonlinear_solution !
00190   _system.get_vector("_old_nonlinear_solution").localize
00191     (*old_local_nonlinear_solution,
00192      _system.get_dof_map().get_send_list());
00193 }

void libMesh::UnsteadySolver::advance_timestep (  )  [virtual, inherited]

This method advances the solution to the next timestep, after a solve() has been performed. Often this will be done after every UnsteadySolver::solve(), but adaptive mesh refinement and/or adaptive time step selection may require some solve() steps to be repeated.

Reimplemented from libMesh::TimeSolver.

Reimplemented in libMesh::AdaptiveTimeSolver.

Definition at line 143 of file unsteady_solver.C.

References libMesh::TimeSolver::_system, libMesh::DifferentiableSystem::deltat, libMesh::UnsteadySolver::first_solve, libMesh::System::get_dof_map(), libMesh::DofMap::get_send_list(), libMesh::System::get_vector(), libMesh::NumericVector< T >::localize(), libMesh::UnsteadySolver::old_local_nonlinear_solution, libMesh::System::solution, libMesh::TimeSolver::solution_history, and libMesh::System::time.

Referenced by libMesh::UnsteadySolver::solve().

00144 {
00145   if (!first_solve)
00146     {
00147       // Store the solution, does nothing by default
00148       // User has to attach appropriate solution_history object for this to
00149       // actually store anything anywhere
00150       solution_history->store();
00151 
00152       _system.time += _system.deltat;
00153     }
00154 
00155   NumericVector<Number> &old_nonlinear_soln =
00156   _system.get_vector("_old_nonlinear_solution");
00157   NumericVector<Number> &nonlinear_solution =
00158     *(_system.solution);
00159 
00160   old_nonlinear_soln = nonlinear_solution;
00161 
00162   old_nonlinear_soln.localize
00163     (*old_local_nonlinear_solution,
00164      _system.get_dof_map().get_send_list());
00165 }

virtual void libMesh::TimeSolver::before_timestep (  )  [inline, virtual, inherited]

This method is for subclasses or users to override to do arbitrary processing between timesteps

Definition at line 152 of file time_solver.h.

00152 {}

virtual AutoPtr<DiffSolver>& libMesh::TimeSolver::diff_solver (  )  [inline, virtual, inherited]

An implicit linear or nonlinear solver to use at each timestep.

Reimplemented in libMesh::AdaptiveTimeSolver.

Definition at line 167 of file time_solver.h.

References libMesh::TimeSolver::_diff_solver.

00167 { return _diff_solver; }

void libMesh::ReferenceCounter::disable_print_counter_info (  )  [static, inherited]

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

00107 {
00108   _enable_print_counter = false;
00109   return;
00110 }

Real libMesh::UnsteadySolver::du ( const SystemNorm norm  )  const [virtual, inherited]

Computes the size of ||u^{n+1} - u^{n}|| in some norm.

Note that, while you can always call this function, its result may or may not be very meaningful. For example, if you call this function right after calling advance_timestep() then you'll get a result of zero since old_nonlinear_solution is set equal to nonlinear_solution in this function.

Implements libMesh::TimeSolver.

Definition at line 218 of file unsteady_solver.C.

References libMesh::TimeSolver::_system, libMesh::System::calculate_norm(), libMesh::System::get_vector(), and libMesh::System::solution.

00219 {
00220 
00221   AutoPtr<NumericVector<Number> > solution_copy =
00222     _system.solution->clone();
00223 
00224   solution_copy->add(-1., _system.get_vector("_old_nonlinear_solution"));
00225 
00226   solution_copy->close();
00227 
00228   return _system.calculate_norm(*solution_copy, norm);
00229 }

bool libMesh::EulerSolver::element_residual ( bool  request_jacobian,
DiffContext context 
) [virtual]

This method uses the DifferentiableSystem's element_time_derivative() and element_constraint() to build a full residual on an element. What combination it uses will depend on theta.

Implements libMesh::TimeSolver.

Definition at line 50 of file euler_solver.C.

References libMesh::TimeSolver::_system, libMesh::DenseVector< T >::add(), libMesh::DifferentiableSystem::deltat, libMesh::DiffContext::dof_indices, libMesh::DiffContext::elem_fixed_solution, libMesh::DiffContext::elem_jacobian, libMesh::DiffContext::elem_reinit(), libMesh::DiffContext::elem_residual, libMesh::DiffContext::elem_solution, libMesh::DiffContext::elem_solution_derivative, libMesh::DifferentiablePhysics::element_constraint(), libMesh::DifferentiablePhysics::element_time_derivative(), libMesh::DifferentiablePhysics::eulerian_residual(), libMesh::DiffContext::fixed_solution_derivative, libMesh::DifferentiablePhysics::mass_residual(), libMesh::UnsteadySolver::old_nonlinear_solution(), libMesh::DenseVector< T >::size(), libMesh::DenseMatrix< T >::swap(), libMesh::DenseVector< T >::swap(), theta, libMesh::System::use_fixed_solution, libMesh::DenseVector< T >::zero(), and libMesh::DenseMatrix< T >::zero().

00052 {
00053   unsigned int n_dofs = context.elem_solution.size();
00054 
00055   // Local nonlinear solution at old timestep
00056   DenseVector<Number> old_elem_solution(n_dofs);
00057   for (unsigned int i=0; i != n_dofs; ++i)
00058     old_elem_solution(i) =
00059       old_nonlinear_solution(context.dof_indices[i]);
00060 
00061   // Local nonlinear solution at time t_theta
00062   DenseVector<Number> theta_solution(context.elem_solution);
00063   theta_solution *= theta;
00064   theta_solution.add(1. - theta, old_elem_solution);
00065 
00066   // Technically the elem_solution_derivative is either theta
00067   // or -1.0 in this implementation, but we scale the former part
00068   // ourselves
00069   context.elem_solution_derivative = 1.0;
00070 
00071 // Technically the fixed_solution_derivative is always theta,
00072 // but we're scaling a whole jacobian by theta after these first
00073 // evaluations
00074   context.fixed_solution_derivative = 1.0;
00075 
00076   // If a fixed solution is requested, we'll use theta_solution
00077   if (_system.use_fixed_solution)
00078     context.elem_fixed_solution = theta_solution;
00079 
00080   // Move theta_->elem_, elem_->theta_
00081   context.elem_solution.swap(theta_solution);
00082 
00083   // Move the mesh into place first if necessary
00084   context.elem_reinit(theta);
00085 
00086   // We're going to compute just the change in elem_residual
00087   // (and possibly elem_jacobian), then add back the old values
00088   DenseVector<Number> old_elem_residual(context.elem_residual);
00089   DenseMatrix<Number> old_elem_jacobian;
00090   if (request_jacobian)
00091     {
00092       old_elem_jacobian = context.elem_jacobian;
00093       context.elem_jacobian.zero();
00094     }
00095   context.elem_residual.zero();
00096 
00097   // Get the time derivative at t_theta
00098   bool jacobian_computed =
00099     _system.element_time_derivative(request_jacobian, context);
00100 
00101   // For a moving mesh problem we may need the pseudoconvection term too
00102   jacobian_computed =
00103     _system.eulerian_residual(jacobian_computed, context) && jacobian_computed;
00104 
00105   // Scale the time-dependent residual and jacobian correctly
00106   context.elem_residual *= _system.deltat;
00107   if (jacobian_computed)
00108     context.elem_jacobian *= (theta * _system.deltat);
00109 
00110   // The fixed_solution_derivative is always theta,
00111   // and now we're done scaling jacobians
00112   context.fixed_solution_derivative = theta;
00113 
00114   // We evaluate mass_residual with the change in solution
00115   // to get the mass matrix, reusing old_elem_solution to hold that
00116   // delta_solution.  We're solving dt*F(u) - du = 0, so our
00117   // delta_solution is old_solution - new_solution.
00118   // We're still keeping elem_solution in theta_solution for now
00119   old_elem_solution -= theta_solution;
00120 
00121   // Move old_->elem_, theta_->old_
00122   context.elem_solution.swap(old_elem_solution);
00123 
00124   // We do a trick here to avoid using a non-1
00125   // elem_solution_derivative:
00126   context.elem_jacobian *= -1.0;
00127   jacobian_computed = _system.mass_residual(jacobian_computed, context) &&
00128     jacobian_computed;
00129   context.elem_jacobian *= -1.0;
00130 
00131   // Move elem_->elem_, old_->theta_
00132   context.elem_solution.swap(theta_solution);
00133 
00134   // Restore the elem position if necessary
00135   context.elem_reinit(1.);
00136 
00137   // Add the constraint term
00138   jacobian_computed = _system.element_constraint(jacobian_computed, context) &&
00139     jacobian_computed;
00140 
00141   // Add back the old residual and jacobian
00142   context.elem_residual += old_elem_residual;
00143   if (request_jacobian)
00144     {
00145       if (jacobian_computed)
00146         context.elem_jacobian += old_elem_jacobian;
00147       else
00148         context.elem_jacobian.swap(old_elem_jacobian);
00149     }
00150 
00151   return jacobian_computed;
00152 }

void libMesh::ReferenceCounter::enable_print_counter_info (  )  [static, inherited]

Methods to enable/disable the reference counter output from print_info()

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

00101 {
00102   _enable_print_counter = true;
00103   return;
00104 }

Real libMesh::EulerSolver::error_order (  )  const [virtual]

Error convergence order: 2 for Crank-Nicolson, 1 otherwise

Implements libMesh::UnsteadySolver.

Definition at line 40 of file euler_solver.C.

References theta.

00041 {
00042   if (theta == 0.5)
00043     return 2.;
00044   return 1.;
00045 }

std::string libMesh::ReferenceCounter::get_info (  )  [static, inherited]

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

References libMesh::ReferenceCounter::_counts, and libMesh::Quality::name().

Referenced by libMesh::ReferenceCounter::print_info().

00048 {
00049 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
00050 
00051   std::ostringstream oss;
00052 
00053   oss << '\n'
00054       << " ---------------------------------------------------------------------------- \n"
00055       << "| Reference count information                                                |\n"
00056       << " ---------------------------------------------------------------------------- \n";
00057 
00058   for (Counts::iterator it = _counts.begin();
00059        it != _counts.end(); ++it)
00060     {
00061       const std::string name(it->first);
00062       const unsigned int creations    = it->second.first;
00063       const unsigned int destructions = it->second.second;
00064 
00065       oss << "| " << name << " reference count information:\n"
00066           << "|  Creations:    " << creations    << '\n'
00067           << "|  Destructions: " << destructions << '\n';
00068     }
00069 
00070   oss << " ---------------------------------------------------------------------------- \n";
00071 
00072   return oss.str();
00073 
00074 #else
00075 
00076   return "";
00077 
00078 #endif
00079 }

void libMesh::ReferenceCounter::increment_constructor_count ( const std::string &  name  )  [inline, protected, inherited]

Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.

Definition at line 163 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

00164 {
00165   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00166   std::pair<unsigned int, unsigned int>& p = _counts[name];
00167 
00168   p.first++;
00169 }

void libMesh::ReferenceCounter::increment_destructor_count ( const std::string &  name  )  [inline, protected, inherited]

Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.

Definition at line 176 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

00177 {
00178   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00179   std::pair<unsigned int, unsigned int>& p = _counts[name];
00180 
00181   p.second++;
00182 }

void libMesh::UnsteadySolver::init (  )  [virtual, inherited]

The initialization function. This method is used to initialize internal data structures before a simulation begins.

Reimplemented from libMesh::TimeSolver.

Reimplemented in libMesh::AdaptiveTimeSolver.

Definition at line 46 of file unsteady_solver.C.

References libMesh::TimeSolver::_system, and libMesh::System::add_vector().

00047 {
00048   TimeSolver::init();
00049 
00050   _system.add_vector("_old_nonlinear_solution");
00051 }

void libMesh::UnsteadySolver::init_data (  )  [virtual, inherited]

The data initialization function. This method is used to initialize internal data structures after the underlying System has been initialized

Reimplemented from libMesh::TimeSolver.

Definition at line 55 of file unsteady_solver.C.

References libMesh::TimeSolver::_system, libMesh::System::get_dof_map(), libMesh::DofMap::get_send_list(), libMeshEnums::GHOSTED, libMesh::System::n_dofs(), libMesh::System::n_local_dofs(), libMesh::UnsteadySolver::old_local_nonlinear_solution, and libMeshEnums::SERIAL.

00056 {
00057 #ifdef LIBMESH_ENABLE_GHOSTED
00058   old_local_nonlinear_solution->init (_system.n_dofs(), _system.n_local_dofs(),
00059                                       _system.get_dof_map().get_send_list(), false,
00060                                       GHOSTED);
00061 #else
00062   old_local_nonlinear_solution->init (_system.n_dofs(), false, SERIAL);
00063 #endif
00064 }

bool libMesh::TimeSolver::is_adjoint (  )  const [inline, inherited]

Accessor for querying whether we need to do a primal or adjoint solve

Definition at line 217 of file time_solver.h.

References libMesh::TimeSolver::_is_adjoint.

Referenced by libMesh::FEMSystem::build_context().

00218   { return _is_adjoint; }

virtual bool libMesh::UnsteadySolver::is_steady (  )  const [inline, virtual, inherited]

This is not a steady-state solver.

Implements libMesh::TimeSolver.

Definition at line 149 of file unsteady_solver.h.

00149 { return false; }

virtual AutoPtr<LinearSolver<Number> >& libMesh::TimeSolver::linear_solver (  )  [inline, virtual, inherited]

An implicit linear solver to use for adjoint and sensitivity problems.

Definition at line 172 of file time_solver.h.

References libMesh::TimeSolver::_linear_solver.

00172 { return _linear_solver; }

static unsigned int libMesh::ReferenceCounter::n_objects (  )  [inline, static, inherited]

Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 79 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

00080   { return _n_objects; }

Number libMesh::UnsteadySolver::old_nonlinear_solution ( const dof_id_type  global_dof_number  )  const [inherited]
Returns:
the old nonlinear solution for the specified global DOF.

Definition at line 207 of file unsteady_solver.C.

References libMesh::TimeSolver::_system, libMesh::System::get_dof_map(), libMesh::DofMap::n_dofs(), and libMesh::UnsteadySolver::old_local_nonlinear_solution.

Referenced by element_residual(), libMesh::Euler2Solver::element_residual(), side_residual(), and libMesh::Euler2Solver::side_residual().

00209 {
00210   libmesh_assert_less (global_dof_number, _system.get_dof_map().n_dofs());
00211   libmesh_assert_less (global_dof_number, old_local_nonlinear_solution->size());
00212 
00213   return (*old_local_nonlinear_solution)(global_dof_number);
00214 }

void libMesh::ReferenceCounter::print_info ( std::ostream &  out = libMesh::out  )  [static, inherited]

Prints the reference information, by default to libMesh::out.

Definition at line 88 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::get_info().

00089 {
00090   if( _enable_print_counter ) out_stream << ReferenceCounter::get_info();
00091 }

void libMesh::UnsteadySolver::reinit (  )  [virtual, inherited]

The reinitialization function. This method is used to resize internal data vectors after a mesh change.

Reimplemented from libMesh::TimeSolver.

Reimplemented in libMesh::AdaptiveTimeSolver.

Definition at line 68 of file unsteady_solver.C.

References libMesh::TimeSolver::_system, libMesh::System::get_dof_map(), libMesh::DofMap::get_send_list(), libMeshEnums::GHOSTED, libMesh::System::n_dofs(), libMesh::System::n_local_dofs(), libMesh::UnsteadySolver::old_local_nonlinear_solution, and libMeshEnums::SERIAL.

00069 {
00070   TimeSolver::reinit();
00071 
00072 #ifdef LIBMESH_ENABLE_GHOSTED
00073   old_local_nonlinear_solution->init (_system.n_dofs(), _system.n_local_dofs(),
00074                                       _system.get_dof_map().get_send_list(), false,
00075                                       GHOSTED);
00076 #else
00077   old_local_nonlinear_solution->init (_system.n_dofs(), false, SERIAL);
00078 #endif
00079 
00080 }

void libMesh::UnsteadySolver::retrieve_timestep (  )  [virtual, inherited]

This method retrieves all the stored solutions at the current system.time

Reimplemented from libMesh::TimeSolver.

Definition at line 195 of file unsteady_solver.C.

References libMesh::TimeSolver::_system, libMesh::System::get_dof_map(), libMesh::DofMap::get_send_list(), libMesh::System::get_vector(), libMesh::NumericVector< T >::localize(), libMesh::UnsteadySolver::old_local_nonlinear_solution, and libMesh::TimeSolver::solution_history.

00196   {
00197     // Retrieve all the stored vectors at the current time
00198     solution_history->retrieve();
00199 
00200     // Dont forget to localize the old_nonlinear_solution !
00201     _system.get_vector("_old_nonlinear_solution").localize
00202     (*old_local_nonlinear_solution,
00203      _system.get_dof_map().get_send_list());
00204   }

void libMesh::TimeSolver::set_is_adjoint ( bool  _is_adjoint_value  )  [inline, inherited]

Accessor for setting whether we need to do a primal or adjoint solve

Definition at line 224 of file time_solver.h.

References libMesh::TimeSolver::_is_adjoint.

Referenced by libMesh::DifferentiableSystem::adjoint_solve(), libMesh::FEMSystem::postprocess(), and libMesh::DifferentiableSystem::solve().

00225   { _is_adjoint = _is_adjoint_value; }

void libMesh::TimeSolver::set_solution_history ( const SolutionHistory _solution_history  )  [inherited]

A setter function users will employ if they need to do something other than save no solution history

Definition at line 91 of file time_solver.C.

References libMesh::SolutionHistory::clone(), and libMesh::TimeSolver::solution_history.

00092  {
00093    solution_history = _solution_history.clone();
00094  }

bool libMesh::EulerSolver::side_residual ( bool  request_jacobian,
DiffContext context 
) [virtual]

This method uses the DifferentiableSystem's side_time_derivative() and side_constraint() to build a full residual on an element's side. What combination it uses will depend on theta.

Implements libMesh::TimeSolver.

Definition at line 156 of file euler_solver.C.

References libMesh::TimeSolver::_system, libMesh::DenseVector< T >::add(), libMesh::DifferentiableSystem::deltat, libMesh::DiffContext::dof_indices, libMesh::DiffContext::elem_fixed_solution, libMesh::DiffContext::elem_jacobian, libMesh::DiffContext::elem_residual, libMesh::DiffContext::elem_side_reinit(), libMesh::DiffContext::elem_solution, libMesh::DiffContext::elem_solution_derivative, libMesh::DiffContext::fixed_solution_derivative, libMesh::UnsteadySolver::old_nonlinear_solution(), libMesh::DifferentiablePhysics::side_constraint(), libMesh::DifferentiablePhysics::side_mass_residual(), libMesh::DifferentiablePhysics::side_time_derivative(), libMesh::DenseVector< T >::size(), libMesh::DenseMatrix< T >::swap(), libMesh::DenseVector< T >::swap(), theta, libMesh::System::use_fixed_solution, libMesh::DenseVector< T >::zero(), and libMesh::DenseMatrix< T >::zero().

00158 {
00159   unsigned int n_dofs = context.elem_solution.size();
00160 
00161   // Local nonlinear solution at old timestep
00162   DenseVector<Number> old_elem_solution(n_dofs);
00163   for (unsigned int i=0; i != n_dofs; ++i)
00164     old_elem_solution(i) =
00165       old_nonlinear_solution(context.dof_indices[i]);
00166 
00167   // Local nonlinear solution at time t_theta
00168   DenseVector<Number> theta_solution(context.elem_solution);
00169   theta_solution *= theta;
00170   theta_solution.add(1. - theta, old_elem_solution);
00171 
00172   // Technically the elem_solution_derivative is either theta
00173   // or 1.0 in this implementation, but we scale the former part
00174   // ourselves
00175   context.elem_solution_derivative = 1.0;
00176 
00177 // Technically the fixed_solution_derivative is always theta,
00178 // but we're scaling a whole jacobian by theta after these first
00179 // evaluations
00180   context.fixed_solution_derivative = 1.0;
00181 
00182   // If a fixed solution is requested, we'll use theta_solution
00183   if (_system.use_fixed_solution)
00184     context.elem_fixed_solution = theta_solution;
00185 
00186   // Move theta_->elem_, elem_->theta_
00187   context.elem_solution.swap(theta_solution);
00188 
00189   // Move the mesh into place first if necessary
00190   context.elem_side_reinit(theta);
00191 
00192   // We're going to compute just the change in elem_residual
00193   // (and possibly elem_jacobian), then add back the old values
00194   DenseVector<Number> old_elem_residual(context.elem_residual);
00195   DenseMatrix<Number> old_elem_jacobian;
00196   if (request_jacobian)
00197     {
00198       old_elem_jacobian = context.elem_jacobian;
00199       context.elem_jacobian.zero();
00200     }
00201   context.elem_residual.zero();
00202 
00203   // Get the time derivative at t_theta
00204   bool jacobian_computed =
00205     _system.side_time_derivative(request_jacobian, context);
00206 
00207   // Scale the time-dependent residual and jacobian correctly
00208   context.elem_residual *= _system.deltat;
00209   if (jacobian_computed)
00210     context.elem_jacobian *= (theta * _system.deltat);
00211 
00212   // The fixed_solution_derivative is always theta,
00213   // and now we're done scaling jacobians
00214   context.fixed_solution_derivative = theta;
00215 
00216   // We evaluate side_mass_residual with the change in solution
00217   // to get the mass matrix, reusing old_elem_solution to hold that
00218   // delta_solution.  We're solving dt*F(u) - du = 0, so our
00219   // delta_solution is old_solution - new_solution.
00220   // We're still keeping elem_solution in theta_solution for now
00221   old_elem_solution -= theta_solution;
00222 
00223   // Move old_->elem_, theta_->old_
00224   context.elem_solution.swap(old_elem_solution);
00225 
00226   // We do a trick here to avoid using a non-1
00227   // elem_solution_derivative:
00228   context.elem_jacobian *= -1.0;
00229   jacobian_computed = _system.side_mass_residual(jacobian_computed, context) &&
00230     jacobian_computed;
00231   context.elem_jacobian *= -1.0;
00232 
00233   // Move elem_->elem_, old_->theta_
00234   context.elem_solution.swap(theta_solution);
00235 
00236   // Restore the elem position if necessary
00237   context.elem_side_reinit(1.);
00238 
00239   // Add the constraint term
00240   jacobian_computed = _system.side_constraint(jacobian_computed, context) &&
00241     jacobian_computed;
00242 
00243   // Add back the old residual and jacobian
00244   context.elem_residual += old_elem_residual;
00245   if (request_jacobian)
00246     {
00247       if (jacobian_computed)
00248         context.elem_jacobian += old_elem_jacobian;
00249       else
00250         context.elem_jacobian.swap(old_elem_jacobian);
00251     }
00252 
00253   return jacobian_computed;
00254 }

void libMesh::UnsteadySolver::solve (  )  [virtual, inherited]

This method solves for the solution at the next timestep. Usually we will only need to solve one (non)linear system per timestep, but more complex subclasses may override this.

Reimplemented from libMesh::TimeSolver.

Reimplemented in libMesh::AdaptiveTimeSolver, and libMesh::TwostepTimeSolver.

Definition at line 84 of file unsteady_solver.C.

References libMesh::TimeSolver::_diff_solver, libMesh::TimeSolver::_system, libMesh::UnsteadySolver::advance_timestep(), libMesh::DifferentiableSystem::deltat, libMesh::DiffSolver::DIVERGED_BACKTRACKING_FAILURE, libMesh::DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS, libMesh::UnsteadySolver::first_solve, libMesh::out, libMesh::TimeSolver::quiet, and libMesh::TimeSolver::reduce_deltat_on_diffsolver_failure.

00085 {
00086   if (first_solve)
00087     {
00088       advance_timestep();
00089       first_solve = false;
00090     }
00091 
00092   unsigned int solve_result = _diff_solver->solve();
00093 
00094   // If we requested the UnsteadySolver to attempt reducing dt after a
00095   // failed DiffSolver solve, check the results of the solve now.
00096   if (reduce_deltat_on_diffsolver_failure)
00097     {
00098       bool backtracking_failed =
00099         solve_result & DiffSolver::DIVERGED_BACKTRACKING_FAILURE;
00100 
00101       bool max_iterations =
00102         solve_result & DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS;
00103 
00104       if (backtracking_failed || max_iterations)
00105         {
00106           // Cut timestep in half
00107           for (unsigned int nr=0; nr<reduce_deltat_on_diffsolver_failure; ++nr)
00108             {
00109               _system.deltat *= 0.5;
00110               libMesh::out << "Newton backtracking failed.  Trying with smaller timestep, dt="
00111                             << _system.deltat << std::endl;
00112 
00113               solve_result = _diff_solver->solve();
00114 
00115               // Check solve results with reduced timestep
00116               bool backtracking_still_failed =
00117                 solve_result & DiffSolver::DIVERGED_BACKTRACKING_FAILURE;
00118 
00119               bool backtracking_max_iterations =
00120                 solve_result & DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS;
00121 
00122               if (!backtracking_still_failed && !backtracking_max_iterations)
00123                 {
00124                   if (!quiet)
00125                     libMesh::out << "Reduced dt solve succeeded." << std::endl;
00126                   return;
00127                 }
00128             }
00129 
00130           // If we made it here, we still couldn't converge the solve after
00131           // reducing deltat
00132           libMesh::out << "DiffSolver::solve() did not succeed after "
00133                         << reduce_deltat_on_diffsolver_failure
00134                         << " attempts." << std::endl;
00135           libmesh_convergence_failure();
00136 
00137         } // end if (backtracking_failed || max_iterations)
00138     } // end if (reduce_deltat_on_diffsolver_failure)
00139 }

sys_type& libMesh::TimeSolver::system (  )  [inline, inherited]
Returns:
a writeable reference to the system we are solving.

Definition at line 162 of file time_solver.h.

References libMesh::TimeSolver::_system.

00162 { return _system; }

const sys_type& libMesh::TimeSolver::system (  )  const [inline, inherited]
Returns:
a constant reference to the system we are solving.

Definition at line 157 of file time_solver.h.

References libMesh::TimeSolver::_system.

Referenced by libMesh::TimeSolver::reinit(), and libMesh::TimeSolver::solve().

00157 { return _system; }


Member Data Documentation

An implicit linear or nonlinear solver to use at each timestep.

Definition at line 232 of file time_solver.h.

Referenced by libMesh::TimeSolver::diff_solver(), libMesh::TimeSolver::init(), libMesh::TimeSolver::reinit(), libMesh::UnsteadySolver::solve(), and libMesh::TimeSolver::solve().

bool libMesh::ReferenceCounter::_enable_print_counter = true [static, protected, inherited]

Flag to control whether reference count information is printed when print_info is called.

Definition at line 137 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::disable_print_counter_info(), libMesh::ReferenceCounter::enable_print_counter_info(), and libMesh::ReferenceCounter::print_info().

An implicit linear solver to use for adjoint problems.

Definition at line 237 of file time_solver.h.

Referenced by libMesh::TimeSolver::init(), libMesh::TimeSolver::linear_solver(), and libMesh::TimeSolver::reinit().

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 131 of file reference_counter.h.

Threads::atomic< unsigned int > libMesh::ReferenceCounter::_n_objects [static, protected, inherited]

The number of objects. Print the reference count information when the number returns to 0.

Definition at line 126 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::n_objects(), libMesh::ReferenceCounter::ReferenceCounter(), and libMesh::ReferenceCounter::~ReferenceCounter().

bool libMesh::UnsteadySolver::first_adjoint_step [protected, inherited]

A bool that will be true the first time adjoint_advance_timestep() is called, (when the primal solution is to be used to set adjoint boundary conditions) and false thereafter

Definition at line 163 of file unsteady_solver.h.

Referenced by libMesh::UnsteadySolver::adjoint_advance_timestep().

bool libMesh::UnsteadySolver::first_solve [protected, inherited]

A bool that will be true the first time solve() is called, and false thereafter

Reimplemented from libMesh::TimeSolver.

Definition at line 157 of file unsteady_solver.h.

Referenced by libMesh::UnsteadySolver::advance_timestep(), libMesh::AdaptiveTimeSolver::advance_timestep(), libMesh::UnsteadySolver::solve(), and libMesh::TwostepTimeSolver::solve().

bool libMesh::TimeSolver::quiet [inherited]

Print extra debugging information if quiet == false.

Definition at line 177 of file time_solver.h.

Referenced by libMesh::UnsteadySolver::solve(), libMesh::TwostepTimeSolver::solve(), and libMesh::EigenTimeSolver::solve().

This value (which defaults to zero) is the number of times the TimeSolver is allowed to halve deltat and let the DiffSolver repeat the latest failed solve with a reduced timestep. Note that this has no effect for SteadySolvers. Note that you must set at least one of the DiffSolver flags "continue_after_max_iterations" or "continue_after_backtrack_failure" to allow the TimeSolver to retry the solve.

Definition at line 205 of file time_solver.h.

Referenced by libMesh::UnsteadySolver::solve(), and libMesh::TwostepTimeSolver::solve().

The value for the theta method to employ: 1.0 corresponds to backwards Euler, 0.0 corresponds to forwards Euler, 0.5 corresponds to Crank-Nicolson.

Definition at line 93 of file euler_solver.h.

Referenced by element_residual(), error_order(), and side_residual().


The documentation for this class was generated from the following files:

Site Created By: libMesh Developers
Last modified: February 05 2013 19:55:18 UTC

Hosted By:
SourceForge.net Logo