/[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 4818 by caltinay, Mon Mar 31 00:16:20 2014 UTC revision 4836 by caltinay, Mon Apr 7 05:51:55 2014 UTC
# Line 92  void Paso_FCT_Solver_free(Paso_FCT_Solve Line 92  void Paso_FCT_Solver_free(Paso_FCT_Solve
92    
93  double Paso_FCT_Solver_getSafeTimeStepSize(Paso_TransportProblem* fctp)  double Paso_FCT_Solver_getSafeTimeStepSize(Paso_TransportProblem* fctp)
94  {  {
95     dim_t i, n;     dim_t i;
96     double dt_max_loc=LARGE_POSITIVE_FLOAT;     double dt_max_loc=LARGE_POSITIVE_FLOAT;
97     double dt_max=LARGE_POSITIVE_FLOAT;     double dt_max=LARGE_POSITIVE_FLOAT;
98     n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     const dim_t n = fctp->transport_matrix->getTotalNumRows();
99     /* set low order transport operator */     /* set low order transport operator */
100     Paso_FCT_setLowOrderOperator(fctp);     Paso_FCT_setLowOrderOperator(fctp);
101                        
# Line 136  void Paso_FCT_Solver_initialize(const do Line 136  void Paso_FCT_Solver_initialize(const do
136  {  {
137     Paso_TransportProblem* fctp = fct_solver->transportproblem;     Paso_TransportProblem* fctp = fct_solver->transportproblem;
138     const index_t* main_iptr=Paso_TransportProblem_borrowMainDiagonalPointer(fctp);     const index_t* main_iptr=Paso_TransportProblem_borrowMainDiagonalPointer(fctp);
139     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     const dim_t n = fctp->transport_matrix->getTotalNumRows();
140     const double theta = Paso_FCT_Solver_getTheta(fct_solver);     const double theta = Paso_FCT_Solver_getTheta(fct_solver);
141     const double omega=1./(dt* theta);     const double omega=1./(dt* theta);
142     dim_t i;     dim_t i;
# Line 144  void Paso_FCT_Solver_initialize(const do Line 144  void Paso_FCT_Solver_initialize(const do
144        
145    
146    
147     Paso_solve_free(fctp->iteration_matrix);     Paso_solve_free(fctp->iteration_matrix.get());
148     /*     /*
149      *   fctp->iteration_matrix[i,i]=m[i]/(dt theta) -l[i,i]      *   fctp->iteration_matrix[i,i]=m[i]/(dt theta) -l[i,i]
150      *      *
# Line 175  void Paso_FCT_Solver_initialize(const do Line 175  void Paso_FCT_Solver_initialize(const do
175      options2.sweeps=-1;      options2.sweeps=-1;
176            
177      Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);      Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
178      Paso_SystemMatrix_setPreconditioner(fctp->iteration_matrix, &options2);      fctp->iteration_matrix->setPreconditioner(&options2);
179      Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);      Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
180  }  }
181    
# Line 198  err_t Paso_FCT_Solver_update(Paso_FCT_So Line 198  err_t Paso_FCT_Solver_update(Paso_FCT_So
198          err_out = SOLVER_INPUT_ERROR;          err_out = SOLVER_INPUT_ERROR;
199      }      }
200      return err_out;      return err_out;
     
201  }  }
202    
203  /* linear crank-nicolson update */  /* linear crank-nicolson update */
# Line 209  err_t Paso_FCT_Solver_update_LCN(Paso_FC Line 208  err_t Paso_FCT_Solver_update_LCN(Paso_FC
208      double *b = fct_solver->b;      double *b = fct_solver->b;
209      double const RTOL = options->tolerance;      double const RTOL = options->tolerance;
210      const dim_t n=Paso_TransportProblem_getTotalNumRows(fct_solver->transportproblem);      const dim_t n=Paso_TransportProblem_getTotalNumRows(fct_solver->transportproblem);
211      Paso_SystemMatrix * iteration_matrix = fct_solver->transportproblem->iteration_matrix;      paso::SystemMatrix_ptr iteration_matrix(fct_solver->transportproblem->iteration_matrix);
212      const index_t*  main_iptr=Paso_TransportProblem_borrowMainDiagonalPointer(fct_solver->transportproblem);      const index_t*  main_iptr=Paso_TransportProblem_borrowMainDiagonalPointer(fct_solver->transportproblem);
213      err_t errorCode = SOLVER_NO_ERROR;      err_t errorCode = SOLVER_NO_ERROR;
214      double norm_u_tilde;      double norm_u_tilde;
# Line 283  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 282  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
282     dim_t i;     dim_t i;
283     double norm_u_tilde, ATOL, norm_du=LARGE_POSITIVE_FLOAT, norm_du_old, rate=1.;     double norm_u_tilde, ATOL, norm_du=LARGE_POSITIVE_FLOAT, norm_du_old, rate=1.;
284     err_t errorCode=SOLVER_NO_ERROR;     err_t errorCode=SOLVER_NO_ERROR;
285     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     const dim_t n = fctp->transport_matrix->getTotalNumRows();
286     const double atol=options->absolute_tolerance;       const double atol=options->absolute_tolerance;  
287     const double rtol=options->tolerance;     const double rtol=options->tolerance;
288     const dim_t max_m=options->iter_max;     const dim_t max_m=options->iter_max;
# Line 437  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 436  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
436   *   for BE : theta = 1.   *   for BE : theta = 1.
437   */   */
438    
439  void Paso_FCT_setAntiDiffusionFlux_CN(Paso_SystemMatrix *flux_matrix,  void Paso_FCT_setAntiDiffusionFlux_CN(paso::SystemMatrix_ptr flux_matrix,
440                                        const Paso_TransportProblem* fct,                                        const Paso_TransportProblem* fct,
441                                        const double dt,                                        const double dt,
442                                        paso::const_Coupler_ptr u_coupler,                                          paso::const_Coupler_ptr u_coupler,  
# Line 452  void Paso_FCT_setAntiDiffusionFlux_CN(Pa Line 451  void Paso_FCT_setAntiDiffusionFlux_CN(Pa
451    const double *remote_u_old=u_old_coupler->borrowRemoteData();    const double *remote_u_old=u_old_coupler->borrowRemoteData();
452    const double dt_half= dt/2;    const double dt_half= dt/2;
453    paso::const_SystemMatrixPattern_ptr pattern(fct->iteration_matrix->pattern);    paso::const_SystemMatrixPattern_ptr pattern(fct->iteration_matrix->pattern);
454    const dim_t n=Paso_SystemMatrix_getTotalNumRows(fct->iteration_matrix);    const dim_t n = fct->iteration_matrix->getTotalNumRows();
455    
456    #pragma omp parallel for schedule(static) private(i, iptr_ij)    #pragma omp parallel for schedule(static) private(i, iptr_ij)
457    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
# Line 483  void Paso_FCT_setAntiDiffusionFlux_CN(Pa Line 482  void Paso_FCT_setAntiDiffusionFlux_CN(Pa
482    }    }
483  }  }
484    
485  void Paso_FCT_setAntiDiffusionFlux_BE(Paso_SystemMatrix *flux_matrix,  void Paso_FCT_setAntiDiffusionFlux_BE(paso::SystemMatrix_ptr flux_matrix,
486                                        const Paso_TransportProblem* fct,                                        const Paso_TransportProblem* fct,
487                                        const double dt,                                        const double dt,
488                                        paso::const_Coupler_ptr u_coupler,                                          paso::const_Coupler_ptr u_coupler,  
# Line 497  void Paso_FCT_setAntiDiffusionFlux_BE(Pa Line 496  void Paso_FCT_setAntiDiffusionFlux_BE(Pa
496    const double *remote_u = u_coupler->borrowRemoteData();    const double *remote_u = u_coupler->borrowRemoteData();
497    const double *remote_u_old = u_old_coupler->borrowRemoteData();    const double *remote_u_old = u_old_coupler->borrowRemoteData();
498    paso::const_SystemMatrixPattern_ptr pattern(fct->iteration_matrix->pattern);    paso::const_SystemMatrixPattern_ptr pattern(fct->iteration_matrix->pattern);
499    const dim_t  n=Paso_SystemMatrix_getTotalNumRows(fct->iteration_matrix);    const dim_t n = fct->iteration_matrix->getTotalNumRows();
500    
501    #pragma omp parallel for schedule(static) private(i, iptr_ij)    #pragma omp parallel for schedule(static) private(i, iptr_ij)
502    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
# Line 537  void Paso_FCT_setAntiDiffusionFlux_BE(Pa Line 536  void Paso_FCT_setAntiDiffusionFlux_BE(Pa
536   *   *
537   */   */
538        
539  void Paso_FCT_setAntiDiffusionFlux_linearCN(Paso_SystemMatrix *flux_matrix,  void Paso_FCT_setAntiDiffusionFlux_linearCN(paso::SystemMatrix_ptr flux_matrix,
540                          const Paso_TransportProblem* fct, const double dt,                          const Paso_TransportProblem* fct, const double dt,
541                          paso::const_Coupler_ptr u_tilde_coupler,                          paso::const_Coupler_ptr u_tilde_coupler,
542                          paso::const_Coupler_ptr u_old_coupler)                          paso::const_Coupler_ptr u_old_coupler)
# Line 550  void Paso_FCT_setAntiDiffusionFlux_linea Line 549  void Paso_FCT_setAntiDiffusionFlux_linea
549    const double *remote_u_tilde = u_tilde_coupler->borrowRemoteData();    const double *remote_u_tilde = u_tilde_coupler->borrowRemoteData();
550    const double *remote_u_old = u_old_coupler->borrowRemoteData();    const double *remote_u_old = u_old_coupler->borrowRemoteData();
551    paso::const_SystemMatrixPattern_ptr pattern(fct->iteration_matrix->pattern);    paso::const_SystemMatrixPattern_ptr pattern(fct->iteration_matrix->pattern);
552    const dim_t n=Paso_SystemMatrix_getTotalNumRows(fct->iteration_matrix);    const dim_t n = fct->iteration_matrix->getTotalNumRows();
553    
554    #pragma omp parallel for schedule(static) private(i, iptr_ij)    #pragma omp parallel for schedule(static) private(i, iptr_ij)
555    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
# Line 609  void Paso_FCT_setLowOrderOperator(Paso_T Line 608  void Paso_FCT_setLowOrderOperator(Paso_T
608    index_t iptr_ij, iptr_ji;    index_t iptr_ij, iptr_ji;
609    const index_t*  main_iptr=Paso_TransportProblem_borrowMainDiagonalPointer(fc);    const index_t*  main_iptr=Paso_TransportProblem_borrowMainDiagonalPointer(fc);
610    
611    if (fc->iteration_matrix==NULL) {    if (!fc->iteration_matrix.get()) {
612        fc->iteration_matrix=Paso_SystemMatrix_alloc(fc->transport_matrix->type,        fc->iteration_matrix.reset(new paso::SystemMatrix(
613                                                     fc->transport_matrix->pattern,                    fc->transport_matrix->type, fc->transport_matrix->pattern,
614                                                     fc->transport_matrix->row_block_size,                    fc->transport_matrix->row_block_size,
615                                                     fc->transport_matrix->col_block_size, TRUE);                    fc->transport_matrix->col_block_size, true));
616    }    }
617    
618    if (Esys_noError()) {    if (Esys_noError()) {
619        paso::const_SystemMatrixPattern_ptr pattern(fc->iteration_matrix->pattern);        paso::const_SystemMatrixPattern_ptr pattern(fc->iteration_matrix->pattern);
620        const dim_t n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);        const dim_t n = fc->iteration_matrix->getTotalNumRows();
621        #pragma omp parallel for private(i, iptr_ij, iptr_ji)  schedule(static)        #pragma omp parallel for private(i, iptr_ij, iptr_ji)  schedule(static)
622        for (i = 0; i < n; ++i) {        for (i = 0; i < n; ++i) {
623            double sum=fc->transport_matrix->mainBlock->val[main_iptr[i]];            double sum=fc->transport_matrix->mainBlock->val[main_iptr[i]];
# Line 682  void Paso_FCT_Solver_setMuPaLu(double* o Line 681  void Paso_FCT_Solver_setMuPaLu(double* o
681                                 const double* M,                                 const double* M,
682                                 paso::const_Coupler_ptr u_coupler,                                 paso::const_Coupler_ptr u_coupler,
683                                 const double a,                                 const double a,
684                                 const Paso_SystemMatrix *L)                                 paso::const_SystemMatrix_ptr L)
685  {  {
686    dim_t i;    dim_t i;
687    paso::const_SystemMatrixPattern_ptr pattern(L->pattern);    paso::const_SystemMatrixPattern_ptr pattern(L->pattern);
688    const double *u = u_coupler->borrowLocalData();    const double *u = u_coupler->borrowLocalData();
689    const double *remote_u = u_coupler->borrowRemoteData();    const double *remote_u = u_coupler->borrowRemoteData();
690    index_t iptr_ij;    index_t iptr_ij;
691    const dim_t n=Paso_SystemMatrix_getTotalNumRows(L);    const dim_t n = L->getTotalNumRows();
692    
693    #pragma omp parallel for private(i) schedule(static)    #pragma omp parallel for private(i) schedule(static)
694    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {

Legend:
Removed from v.4818  
changed lines
  Added in v.4836

  ViewVC Help
Powered by ViewVC 1.1.26