/[escript]/trunk/paso/src/FCTSolver.c
ViewVC logotype

Diff of /trunk/paso/src/FCTSolver.c

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

revision 3005 by gross, Thu Apr 22 05:59:31 2010 UTC revision 3014 by gross, Wed Apr 28 04:05:21 2010 UTC
# Line 116  double Paso_FCTSolver_getSafeTimeStepSiz Line 116  double Paso_FCTSolver_getSafeTimeStepSiz
116     n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
117    
118    
119       /* set low order transport operator */     /* set low order transport operator */
120       Paso_FCTSolver_setLowOrderOperator(fctp);     Paso_FCTSolver_setLowOrderOperator(fctp);
121                        
122     if (Paso_noError()) {     if (Paso_noError()) {
123          /*          /*
# Line 176  err_t Paso_FCTSolver_Function_call(Paso_ Line 176  err_t Paso_FCTSolver_Function_call(Paso_
176    
177  err_t Paso_FCTSolver_solve(Paso_Function* F, double* u, double dt, Paso_Options* options, Paso_Performance *pp)  err_t Paso_FCTSolver_solve(Paso_Function* F, double* u, double dt, Paso_Options* options, Paso_Performance *pp)
178  {  {
179     double norm_u_tilde, ATOL, norm_u, norm_du, *du=NULL;     const dim_t num_critical_rates_max=3; /* number of rates >=critical_rate accepted before divergence is triggered */
180       const double critical_rate=0.6;   /* expected value of convergence rate */
181      
182       double norm_u_tilde, ATOL, norm_u, norm_du=LARGE_POSITIVE_FLOAT, norm_du_old, *du=NULL, rate;
183     err_t errorCode=SOLVER_NO_ERROR;     err_t errorCode=SOLVER_NO_ERROR;
184     Paso_FCTSolver *more=(Paso_FCTSolver *) (F->more);     Paso_FCTSolver *more=(Paso_FCTSolver *) (F->more);
185     const double omega=1./dt*( more->transportproblem->useBackwardEuler ? 1. : 2.);     const double omega=1./(dt*Paso_Transport_getTheta(more->transportproblem));
186     Paso_TransportProblem* fctp=more->transportproblem;     Paso_TransportProblem* fctp=more->transportproblem;
187     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
188     const double atol=options->absolute_tolerance;       const double atol=options->absolute_tolerance;  
189     const double rtol=options->tolerance;     const double rtol=options->tolerance;
190     const dim_t max_m=options->iter_max;     const dim_t max_m=options->iter_max;
191     dim_t m=0;     dim_t m=0, num_critical_rates=0 ;
192     bool_t converged=FALSE, max_m_reached=FALSE;     bool_t converged=FALSE, max_m_reached=FALSE,diverged=FALSE;
193    
194     /* tolerance? */     /* tolerance? */
195      more->dt=dt;      more->dt=dt;
# Line 205  err_t Paso_FCTSolver_solve(Paso_Function Line 208  err_t Paso_FCTSolver_solve(Paso_Function
208          options->tolerance=rtol;          options->tolerance=rtol;
209      } else {      } else {
210              m=0;              m=0;
             converged=FALSE;  
             max_m_reached=FALSE;  
211          du=MEMALLOC(n,double);          du=MEMALLOC(n,double);
212          if (Paso_checkPtr(du)) {          if (Paso_checkPtr(du)) {
213                     errorCode=SOLVER_MEMORY_ERROR;                     errorCode=SOLVER_MEMORY_ERROR;
214          } else {          } else {
215               /* tolerance? */               /* tolerance? */
216                            
217                   while ( (!converged) && (! max_m_reached) && Paso_noError()) {                   while ( (!converged) && (!diverged) && (! max_m_reached) && Paso_noError()) {
218                  errorCode=Paso_FCTSolver_Function_call(F,du, u, pp);                  errorCode=Paso_FCTSolver_Function_call(F,du, u, pp);
219                      Paso_Update(n,1.,u,omega,du);                      Paso_Update(n,1.,u,omega,du);
220                       norm_u=Paso_lsup(n,u, fctp->mpi_info);                       norm_u=Paso_lsup(n,u, fctp->mpi_info);
221                 norm_du_old=norm_du;
222                       norm_du=Paso_lsup(n,du, fctp->mpi_info);                       norm_du=Paso_lsup(n,du, fctp->mpi_info);
223                       if (options->verbose) printf("Paso_FCTSolver_solve: step %d: increment= %e (tolerance = %e)\n",m+1, norm_du*omega, ATOL);               if (m ==0) {
224                             if (options->verbose) printf("Paso_FCTSolver_solve: step %d: increment= %e\n",m+1, norm_du*omega);
225    
226                 } else {
227                     if (norm_du_old > 0.) {
228                         rate=norm_du/norm_du_old;
229                 } else if (norm_du <= 0.) {
230                     rate=0.;
231                 } else {
232                     rate=LARGE_POSITIVE_FLOAT;
233                 }
234                             if (options->verbose) printf("Paso_FCTSolver_solve: step %d: increment= %e (rate = %e)\n",m+1, norm_du*omega, rate);
235                 num_critical_rates+=( rate<critical_rate ? 0 : 1);
236                 }
237    
238                       max_m_reached=(m>max_m);                       max_m_reached=(m>max_m);
239                 diverged = (num_critical_rates >= num_critical_rates_max);
240                       converged=(norm_du*omega <= ATOL) ;                       converged=(norm_du*omega <= ATOL) ;
241                       m++;                       m++;
242                   }                   }
# Line 229  err_t Paso_FCTSolver_solve(Paso_Function Line 246  err_t Paso_FCTSolver_solve(Paso_Function
246              if (converged) {              if (converged) {
247                  if (options->verbose) printf("Paso_FCTSolver_solve: iteration is completed.\n");                  if (options->verbose) printf("Paso_FCTSolver_solve: iteration is completed.\n");
248                      errorCode=SOLVER_NO_ERROR;                      errorCode=SOLVER_NO_ERROR;
249                  } else if (max_m_reached) {                  } else if (diverged) {
250                if (options->verbose) printf("Paso_FCTSolver_solve: divergence.\n");
251                errorCode=SOLVER_DIVERGENCE;
252            } else if (max_m_reached) {
253                 if (options->verbose) printf("Paso_FCTSolver_solve: maximum number of iteration steps reached.\n");
254              errorCode=SOLVER_MAXITER_REACHED;              errorCode=SOLVER_MAXITER_REACHED;
255          }          }
256          }          }
# Line 246  err_t Paso_FCTSolver_setUpRightHandSide( Line 267  err_t Paso_FCTSolver_setUpRightHandSide(
267                      double *RN_m, Paso_Coupler* RN_m_coupler, double* RP_m, Paso_Coupler* RP_m_coupler,                      double *RN_m, Paso_Coupler* RN_m_coupler, double* RP_m, Paso_Coupler* RP_m_coupler,
268                      Paso_Performance* pp)                      Paso_Performance* pp)
269  {  {
270     const double omega=dt* (fctp->useBackwardEuler ? 1. : 0.5);     const double omega=dt* Paso_Transport_getTheta(fctp);
271     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
272     /* distribute u */     /* distribute u */
273     Paso_Coupler_startCollect(u_m_coupler,u_m);     Paso_Coupler_startCollect(u_m_coupler,u_m);
# Line 290  err_t Paso_FCTSolver_setUpRightHandSide( Line 311  err_t Paso_FCTSolver_setUpRightHandSide(
311     return SOLVER_NO_ERROR;     return SOLVER_NO_ERROR;
312  }  }
313    
314    void Paso_FCTSolver_Function_initialize(const double dt, Paso_TransportProblem* fctp, Paso_Options* options, Paso_Performance* pp)
315    {
316       dim_t i;
317       const index_t* main_iptr=Paso_TransportProblem_borrowMainDiagonalPointer(fctp);
318       const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
319       const double theta = Paso_Transport_getTheta(fctp);
320       const double omega=1./(dt* theta);
321       register double m, rtmp4;
322    
323       /*
324        *   fctp->iteration_matrix[i,i]=m[i]/(dt theta) -l[i,i]
325        *
326        */
327        Paso_solve_free(fctp->iteration_matrix);    
328        #pragma omp parallel for private(i,m,rtmp4)
329        for (i = 0; i < n; ++i) {
330               m=fctp->lumped_mass_matrix[i];
331               rtmp4=m*omega-(fctp->main_diagonal_low_order_transport_matrix[i]);
332               fctp->iteration_matrix->mainBlock->val[main_iptr[i]]=rtmp4;
333        }
334        Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
335        Paso_Solver_setPreconditioner(fctp->iteration_matrix,options);
336        Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
337    }
338    
339  void Paso_FCTSolver_setUp(Paso_TransportProblem* fctp, const double dt, const double *u, double* b, double* uTilde,  void Paso_FCTSolver_setUp(Paso_TransportProblem* fctp, const double dt, const double *u, double* b, double* uTilde,
340                            Paso_Coupler* uTilde_coupler, double *QN, Paso_Coupler* QN_coupler, double *QP, Paso_Coupler* QP_coupler,                            Paso_Coupler* uTilde_coupler, double *QN, Paso_Coupler* QN_coupler, double *QP, Paso_Coupler* QP_coupler,
341                            Paso_Options* options, Paso_Performance* pp)                            Paso_Options* options, Paso_Performance* pp)
342  {  {
343     dim_t i;     dim_t i;
    const index_t* main_iptr=NULL;  
344     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
345     const double theta = (fctp->useBackwardEuler ? 1. : 0.5);     register double m;
    const double omega=1./(dt* theta);  
    register double m, u_tilde_i, rtmp4;  
    main_iptr=Paso_TransportProblem_borrowMainDiagonalPointer(fctp);  
346     /* distribute u */     /* distribute u */
347     Paso_Coupler_startCollect(fctp->u_coupler,u);     Paso_Coupler_startCollect(fctp->u_coupler,u);
348     Paso_Coupler_finishCollect(fctp->u_coupler);     Paso_Coupler_finishCollect(fctp->u_coupler);
# Line 319  void Paso_FCTSolver_setUp(Paso_Transport Line 361  void Paso_FCTSolver_setUp(Paso_Transport
361     /*     /*
362      *   uTilde[i]=b[i]/m[i]      *   uTilde[i]=b[i]/m[i]
363      *      *
364      *   fctp->iteration_matrix[i,i]=m[i]/(dt theta) -l[i,i]      */  
365      *      #pragma omp parallel for private(i,m)
     */  
     Paso_solve_free(fctp->iteration_matrix);      
     #pragma omp parallel for private(i,m,u_tilde_i,rtmp4)  
366      for (i = 0; i < n; ++i) {      for (i = 0; i < n; ++i) {
367             m=fctp->lumped_mass_matrix[i];             m=fctp->lumped_mass_matrix[i];
368             u_tilde_i=b[i]/m;             uTilde[i]=b[i]/m;
369             rtmp4=m*omega-(fctp->main_diagonal_low_order_transport_matrix[i]);      }        
            fctp->iteration_matrix->mainBlock->val[main_iptr[i]]=rtmp4;  
            uTilde[i]=u_tilde_i;  
     }  
     Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);  
     Paso_Solver_setPreconditioner(fctp->iteration_matrix,options);  
     Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);  
           
370      /* distribute uTilde: */      /* distribute uTilde: */
371      Paso_Coupler_startCollect(uTilde_coupler,uTilde);      Paso_Coupler_startCollect(uTilde_coupler,uTilde);
372      Paso_Coupler_finishCollect(uTilde_coupler);      Paso_Coupler_finishCollect(uTilde_coupler);

Legend:
Removed from v.3005  
changed lines
  Added in v.3014

  ViewVC Help
Powered by ViewVC 1.1.26