/[escript]/trunk/paso/src/FCT_Solver.cpp
ViewVC logotype

Diff of /trunk/paso/src/FCT_Solver.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4867 by caltinay, Fri Apr 11 12:29:49 2014 UTC revision 4869 by caltinay, Mon Apr 14 10:39:22 2014 UTC
# Line 76  FCT_Solver::~FCT_Solver() Line 76  FCT_Solver::~FCT_Solver()
76  }  }
77    
78  // modifies the main diagonal of the iteration matrix to introduce new dt  // modifies the main diagonal of the iteration matrix to introduce new dt
79  void FCT_Solver::initialize(double _dt, Options* options, Paso_Performance* pp)  void FCT_Solver::initialize(double _dt, Options* options, Performance* pp)
80  {  {
81      const_TransportProblem_ptr fctp(transportproblem);      const_TransportProblem_ptr fctp(transportproblem);
82      const index_t* main_iptr = fctp->borrowMainDiagonalPointer();      const index_t* main_iptr = fctp->borrowMainDiagonalPointer();
# Line 86  void FCT_Solver::initialize(double _dt, Line 86  void FCT_Solver::initialize(double _dt,
86      dim_t i;      dim_t i;
87      Options options2;      Options options2;
88    
89      Paso_solve_free(fctp->iteration_matrix.get());      solve_free(fctp->iteration_matrix.get());
90      //   fctp->iteration_matrix[i,i]=m[i]/(dt theta) -l[i,i]      //   fctp->iteration_matrix[i,i]=m[i]/(dt theta) -l[i,i]
91      dt = _dt;      dt = _dt;
92      #pragma omp parallel for private(i)      #pragma omp parallel for private(i)
# Line 96  void FCT_Solver::initialize(double _dt, Line 96  void FCT_Solver::initialize(double _dt,
96          if ( m_i > 0 ) {          if ( m_i > 0 ) {
97              fctp->iteration_matrix->mainBlock->val[main_iptr[i]] = m_i * omega - l_ii;              fctp->iteration_matrix->mainBlock->val[main_iptr[i]] = m_i * omega - l_ii;
98          } else {          } else {
99              fctp->iteration_matrix->mainBlock->val[main_iptr[i]] = ABS(m_i * omega - l_ii)/(EPSILON*EPSILON);              fctp->iteration_matrix->mainBlock->val[main_iptr[i]] = std::abs(m_i * omega - l_ii)/(EPSILON*EPSILON);
100          }          }
101      }      }
102    
# Line 111  void FCT_Solver::initialize(double _dt, Line 111  void FCT_Solver::initialize(double _dt,
111      options2.use_local_preconditioner = false;      options2.use_local_preconditioner = false;
112      options2.sweeps = -1;      options2.sweeps = -1;
113    
114      Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);      Performance_startMonitor(pp, PERFORMANCE_PRECONDITIONER_INIT);
115      fctp->iteration_matrix->setPreconditioner(&options2);      fctp->iteration_matrix->setPreconditioner(&options2);
116      Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);      Performance_stopMonitor(pp, PERFORMANCE_PRECONDITIONER_INIT);
117  }  }
118    
119  // entry point for update procedures  // entry point for update procedures
120  err_t FCT_Solver::update(double* u, double* u_old, Options* options, Paso_Performance* pp)  err_t FCT_Solver::update(double* u, double* u_old, Options* options, Performance* pp)
121  {  {
122      err_t err_out = SOLVER_NO_ERROR;      err_t err_out = SOLVER_NO_ERROR;
123    
# Line 134  err_t FCT_Solver::update(double* u, doub Line 134  err_t FCT_Solver::update(double* u, doub
134  }  }
135    
136  /// linear crank-nicolson update  /// linear crank-nicolson update
137  err_t FCT_Solver::updateLCN(double* u, double* u_old, Options* options, Paso_Performance *pp)  err_t FCT_Solver::updateLCN(double* u, double* u_old, Options* options, Performance* pp)
138  {  {
139      dim_t sweep_max, i;      dim_t sweep_max, i;
140      double const RTOL = options->tolerance;      double const RTOL = options->tolerance;
# Line 176  err_t FCT_Solver::updateLCN(double* u, d Line 176  err_t FCT_Solver::updateLCN(double* u, d
176      // initial guess is u<- -u + 2*u_tilde      // initial guess is u<- -u + 2*u_tilde
177      util::update(n, -1., u, 2., flux_limiter->u_tilde);      util::update(n, -1., u, 2., flux_limiter->u_tilde);
178    
179      sweep_max = MAX((int) (- 2 * log(RTOL)/log(2.)-0.5),1);      sweep_max = std::max((int) (- 2 * log(RTOL)/log(2.)-0.5),1);
180      norm_u_tilde = util::lsup(n, flux_limiter->u_tilde, flux_limiter->mpi_info);      norm_u_tilde = util::lsup(n, flux_limiter->u_tilde, flux_limiter->mpi_info);
181      if (options->verbose) {      if (options->verbose) {
182          printf("FCT_Solver::updateLCN: u_tilde lsup = %e (rtol = %e, max. sweeps = %d)\n", norm_u_tilde, RTOL*norm_u_tilde, sweep_max);          printf("FCT_Solver::updateLCN: u_tilde lsup = %e (rtol = %e, max. sweeps = %d)\n", norm_u_tilde, RTOL*norm_u_tilde, sweep_max);
# Line 194  err_t FCT_Solver::updateLCN(double* u, d Line 194  err_t FCT_Solver::updateLCN(double* u, d
194      return errorCode;      return errorCode;
195  }  }
196    
197  err_t FCT_Solver::updateNL(double* u, double *u_old, Options* options,  err_t FCT_Solver::updateNL(double* u, double* u_old, Options* options,
198                             Paso_Performance *pp)                             Performance* pp)
199  {  {
200      // number of rates >=critical_rate accepted before divergence is triggered      // number of rates >=critical_rate accepted before divergence is triggered
201      const dim_t num_critical_rates_max = 3;      const dim_t num_critical_rates_max = 3;
# Line 289  err_t FCT_Solver::updateNL(double* u, do Line 289  err_t FCT_Solver::updateNL(double* u, do
289              if (m==0) {              if (m==0) {
290                  tol *= 0.5;                  tol *= 0.5;
291              } else {              } else {
292                  tol *= MIN(MAX(rate*rate, 1e-2), 0.5);                  tol *= std::min(std::max(rate*rate, 1e-2), 0.5);
293              }              }
294              // use BiCGStab with Jacobi preconditioner ( m - omega * L )              // use BiCGStab with Jacobi preconditioner ( m - omega * L )
295              util::zeroes(n,du);              util::zeroes(n,du);
# Line 545  double FCT_Solver::getSafeTimeStepSize(T Line 545  double FCT_Solver::getSafeTimeStepSize(T
545                  const double m_i = fctp->lumped_mass_matrix[i];                  const double m_i = fctp->lumped_mass_matrix[i];
546                  if (m_i > 0) {                  if (m_i > 0) {
547                      if (l_ii<0)                      if (l_ii<0)
548                          dt_max_loc = MIN(dt_max_loc,m_i/(-l_ii));                          dt_max_loc = std::min(dt_max_loc,m_i/(-l_ii));
549                  }                  }
550              }              }
551              #pragma omp critical              #pragma omp critical
552              {              {
553                  dt_max = MIN(dt_max,dt_max_loc);                  dt_max = std::min(dt_max,dt_max_loc);
554              }              }
555          }          }
556  #ifdef ESYS_MPI  #ifdef ESYS_MPI
# Line 662  void FCT_Solver::setMuPaLu(double* out, Line 662  void FCT_Solver::setMuPaLu(double* out,
662              out[i] = u[i];              out[i] = u[i];
663          }          }
664      }      }
665      if (ABS(a) > 0) {      if (std::abs(a) > 0) {
666  #pragma omp parallel for  #pragma omp parallel for
667          for (dim_t i = 0; i < n; ++i) {          for (dim_t i = 0; i < n; ++i) {
668              if (M[i] > 0.) {              if (M[i] > 0.) {

Legend:
Removed from v.4867  
changed lines
  Added in v.4869

  ViewVC Help
Powered by ViewVC 1.1.26