linear_implicit_timestep.h

Go to the documentation of this file.

/*
 * Copyright (c) 2014-2016:  G-CSC, Goethe University Frankfurt
 * Author: Arne Naegel
 *
 * This file is part of UG4.
 *
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 *
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 *
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 *
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 *
 * This program 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.
 */
 
#ifndef LIMEX_H_
#define LIMEX_H_
 
// extern libraries
#include <vector>
#include <cmath>
 
// ug libraries
#include "common/common.h"
#include "common/util/smart_pointer.h"
 
#include "lib_disc/time_disc/time_disc_interface.h"
#include "lib_disc/assemble_interface.h"
#include "lib_disc/operator/non_linear_operator/assembled_non_linear_operator.h"
#include "lib_disc/operator/linear_operator/assembled_linear_operator.h"
 
#include "lib_algebra/algebra_template_define_helper.h"
#include "lib_algebra/operator/debug_writer.h"
 
namespace ug{
 
 
template <class TAlgebra>
class LinearImplicitEuler
: public ITimeDiscretization<TAlgebra>,
  public DebugWritingObject<TAlgebra>
 
{
public:
  typedef ITimeDiscretization<TAlgebra> base_type;
 
  typedef TAlgebra algebra_type;
 
  typedef typename algebra_type::matrix_type matrix_type;
 
  typedef typename algebra_type::vector_type vector_type;
 
  typedef IDomainDiscretization<algebra_type> domain_discretization_type;
 
public:
  LinearImplicitEuler(SmartPtr<IDomainDiscretization<algebra_type> > spDD)
    : ITimeDiscretization<TAlgebra>(spDD),
      m_pPrevSol(NULL),
      m_dt(0.0),
      m_futureTime(0.0),
      m_spMatrixJDisc(spDD), m_spMatrixJOp(SPNULL), m_bMatrixJNeedsUpdate(true), m_useLinearMode(false),
      m_spGammaDisc(SPNULL), m_spGammaOp(SPNULL), m_bGammaNeedsUpdate(true),
      m_useCachedMatrices(true)
  {}
 
  LinearImplicitEuler(SmartPtr<IDomainDiscretization<algebra_type> > spDefectDisc,
            SmartPtr<IDomainDiscretization<algebra_type> > spMatrixJDisc)
    : ITimeDiscretization<TAlgebra>(spDefectDisc),
      m_pPrevSol(NULL),
      m_dt(0.0),
      m_futureTime(0.0),
      m_spMatrixJDisc(spMatrixJDisc), m_spMatrixJOp(SPNULL), m_bMatrixJNeedsUpdate(true), m_useLinearMode(false), 
      m_spGammaDisc(SPNULL), m_spGammaOp(SPNULL), m_bGammaNeedsUpdate(true),
      m_useCachedMatrices(true)
  {}
 
  LinearImplicitEuler(SmartPtr<IDomainDiscretization<algebra_type> > spDefectDisc,
            SmartPtr<IDomainDiscretization<algebra_type> > spMatrixJDisc,
            SmartPtr<IDomainDiscretization<algebra_type> > spGammaDisc)
    : ITimeDiscretization<TAlgebra>(spDefectDisc),
      m_pPrevSol(NULL),
      m_dt(0.0),
      m_futureTime(0.0),
      m_spMatrixJDisc(spMatrixJDisc), m_spMatrixJOp(SPNULL), m_bMatrixJNeedsUpdate(true), m_useLinearMode(false), 
      m_spGammaDisc(spGammaDisc), m_spGammaOp(SPNULL), m_bGammaNeedsUpdate(true),
      m_useCachedMatrices(true)
  {}
 
  virtual ~LinearImplicitEuler(){};
 
public:
  virtual size_t num_prev_steps() const {return m_prevSteps;}
 
  virtual void prepare_step(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
                            number dt);
 
  virtual void prepare_step_elem(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
                                 number dt, const GridLevel& gl);
 
  virtual void finish_step(SmartPtr<VectorTimeSeries<vector_type> > currSol) {};
 
  virtual void finish_step_elem(SmartPtr<VectorTimeSeries<vector_type> > currSol,
                                const GridLevel& gl);
 
  virtual number future_time() const {return m_futureTime;}
 
public:
 
  void assemble_jacobian(matrix_type& J, const vector_type& u, const GridLevel& gl);
 
  void assemble_defect(vector_type& d, const vector_type& u, const GridLevel& gl);
 
  void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl);
 
  void assemble_rhs(vector_type& b, const vector_type& u, const GridLevel& gl);
 
  void assemble_rhs(vector_type& b, const GridLevel& gl);
 
  void adjust_solution(vector_type& u, const GridLevel& gl);
 
   virtual size_t num_stages() const {return 1;};
   virtual void set_stage(size_t stage) {};
 
 
  void enable_linear_mode() { m_useLinearMode = true; }
  void disable_linear_mode() { m_useLinearMode = false; }
  bool use_linear_mode() const { return m_useLinearMode; }
 
  void enable_matrix_cache() { m_useCachedMatrices = true; }
  void disable_matrix_cache() { m_useCachedMatrices = false; }
  void set_matrix_cache(bool useCache) { m_useCachedMatrices = useCache; }
 
protected:
 
  virtual number update_scaling(std::vector<number>& vSM,
                                    std::vector<number>& vSA,
                                    number dt, number currentTime,
                                    ConstSmartPtr<VectorTimeSeries<vector_type> > prevSol)
 
  {
    //  resize scaling factors
    vSM.resize(1);
    vSM[0] = 1.0;
 
    vSA.resize(1);
    vSA[0] = dt;
 
    return currentTime + dt;
  }
 
  static const size_t m_prevSteps=1;      
  std::vector<number> m_vScaleMass;     
  std::vector<number> m_vScaleStiff;      
 
  SmartPtr<VectorTimeSeries<vector_type> > m_pPrevSol;    
  number m_dt;                
  number m_futureTime;            
 
  // discretization for defect
  using base_type::m_spDomDisc;
 
 
  // constant matrix $$ \tau J$$
  SmartPtr<IDomainDiscretization<algebra_type> > m_spMatrixJDisc;
  SmartPtr<AssembledLinearOperator<algebra_type> > m_spMatrixJOp;   
  bool m_bMatrixJNeedsUpdate;
  bool m_useLinearMode;
 
  // Matrix: $\Gamma[u0, u0']$
  SmartPtr<IDomainDiscretization<algebra_type> > m_spGammaDisc;   
  SmartPtr<AssembledLinearOperator<algebra_type> > m_spGammaOp;     
  bool m_bGammaNeedsUpdate;
 
  SmartPtr<matrix_type> m_spMatrixCacheMk;
 
  bool m_useCachedMatrices;
 
public:
 
  void invalidate()
  {
    m_spMatrixCacheMk = SPNULL;
    m_bMatrixJNeedsUpdate = true;
    invalidate_gamma();
  }
 
  void invalidate_gamma()
  { m_spGammaOp = SPNULL; m_bGammaNeedsUpdate = true; }
 
  void set_gamma_disc(SmartPtr<IDomainDiscretization<algebra_type> > spGammaDisc)
  {m_spGammaDisc = spGammaDisc;}
 
  using DebugWritingObject<TAlgebra>::set_debug;
 
protected:
  using DebugWritingObject<TAlgebra>::write_debug;
  using DebugWritingObject<TAlgebra>::debug_writer;
 
};
 
// end group limex
 
}  // namespace ug
#endif