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

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

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

revision 3822 by gross, Mon Feb 13 06:22:06 2012 UTC revision 3824 by gross, Tue Feb 14 01:22:16 2012 UTC
# Line 97  void Paso_FCT_Solver_free(Paso_FCT_Solve Line 97  void Paso_FCT_Solver_free(Paso_FCT_Solve
97  double Paso_FCT_Solver_getSafeTimeStepSize(Paso_TransportProblem* fctp)  double Paso_FCT_Solver_getSafeTimeStepSize(Paso_TransportProblem* fctp)
98  {  {
99     dim_t i, n;     dim_t i, n;
100       double dt_max_loc=LARGE_POSITIVE_FLOAT;
101     double dt_max=LARGE_POSITIVE_FLOAT;     double dt_max=LARGE_POSITIVE_FLOAT;
102     n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
103     /* set low order transport operator */     /* set low order transport operator */
# Line 107  double Paso_FCT_Solver_getSafeTimeStepSi Line 108  double Paso_FCT_Solver_getSafeTimeStepSi
108           *  calculate time step size:                                                     *  calculate time step size:                                          
109          */          */
110          dt_max=LARGE_POSITIVE_FLOAT;          dt_max=LARGE_POSITIVE_FLOAT;
111          #pragma omp parallel private(i)          #pragma omp parallel private(i, dt_max_loc)
112          {          {
113                 double dt_max_loc=LARGE_POSITIVE_FLOAT;                 dt_max_loc=LARGE_POSITIVE_FLOAT;
114                 #pragma omp for schedule(static)                 #pragma omp for schedule(static)
115                 for (i=0;i<n;++i) {                 for (i=0;i<n;++i) {
116                    const double l_ii=fctp->main_diagonal_low_order_transport_matrix[i];                    const double l_ii=fctp->main_diagonal_low_order_transport_matrix[i];
# Line 208  err_t Paso_FCT_Solver_update(Paso_FCT_So Line 209  err_t Paso_FCT_Solver_update(Paso_FCT_So
209  err_t Paso_FCT_Solver_update_LCN(Paso_FCT_Solver *fct_solver, double * u, double *u_old, Paso_Options* options, Paso_Performance *pp)  err_t Paso_FCT_Solver_update_LCN(Paso_FCT_Solver *fct_solver, double * u, double *u_old, Paso_Options* options, Paso_Performance *pp)
210  {  {
211      double const dt = fct_solver->dt;      double const dt = fct_solver->dt;
212      dim_t sweep_max;      dim_t sweep_max, i;
213      double *b = fct_solver->b;      double *b = fct_solver->b;
214      double const RTOL = options->tolerance;      double const RTOL = options->tolerance;
215      const dim_t n=Paso_TransportProblem_getTotalNumRows(fct_solver->transportproblem);      const dim_t n=Paso_TransportProblem_getTotalNumRows(fct_solver->transportproblem);
216      Paso_SystemMatrix * iteration_matrix = fct_solver->transportproblem->iteration_matrix;      Paso_SystemMatrix * iteration_matrix = fct_solver->transportproblem->iteration_matrix;
217        const index_t*  main_iptr=Paso_TransportProblem_borrowMainDiagonalPointer(fct_solver->transportproblem);
218      err_t errorCode = SOLVER_NO_ERROR;      err_t errorCode = SOLVER_NO_ERROR;
219        double norm_u_tilde;
220        
221      Paso_Coupler_startCollect(fct_solver->u_old_coupler,u_old);      Paso_Coupler_startCollect(fct_solver->u_old_coupler,u_old);
222      Paso_Coupler_finishCollect(fct_solver->u_old_coupler);      Paso_Coupler_finishCollect(fct_solver->u_old_coupler);
# Line 240  err_t Paso_FCT_Solver_update_LCN(Paso_FC Line 243  err_t Paso_FCT_Solver_update_LCN(Paso_FC
243     Paso_FCT_FluxLimiter_addLimitedFluxes_Start(fct_solver->flux_limiter);     Paso_FCT_FluxLimiter_addLimitedFluxes_Start(fct_solver->flux_limiter);
244     Paso_FCT_FluxLimiter_addLimitedFluxes_Complete(fct_solver->flux_limiter, b);     Paso_FCT_FluxLimiter_addLimitedFluxes_Complete(fct_solver->flux_limiter, b);
245        
246         Paso_Scale(n, b,fct_solver->omega );    
247     /* solve (m-dt/2*L) u = b in the form (omega*m-L) u = b * omega with omega*dt/2=1 */     /* solve (m-dt/2*L) u = b in the form (omega*m-L) u = b * omega with omega*dt/2=1 */
248         #pragma omp for private(i) schedule(static)
249       for (i = 0; i < n; ++i) {
250           if (! (fct_solver->transportproblem->lumped_mass_matrix[i] > 0 ) ) {
251         b[i] = fct_solver->flux_limiter->u_tilde[i]
252             * fct_solver->transportproblem->iteration_matrix->mainBlock->val[main_iptr[i]];
253           }
254       }
255     /* initial guess is u<- -u + 2*u_tilde */     /* initial guess is u<- -u + 2*u_tilde */
256     Paso_Update(n, -1., u, 2., fct_solver->flux_limiter->u_tilde);     Paso_Update(n, -1., u, 2., fct_solver->flux_limiter->u_tilde);
    Paso_Scale(n, b,fct_solver->omega );  
    sweep_max = MAX((int) (- 2 * log(RTOL)/log(2.)-0.5),1);  
257    
258       sweep_max = MAX((int) (- 2 * log(RTOL)/log(2.)-0.5),1);
259       norm_u_tilde=Paso_lsup(n, fct_solver->flux_limiter->u_tilde, fct_solver->flux_limiter->mpi_info);
260     if (options->verbose) {     if (options->verbose) {
261         const double norm_u_tilde=Paso_lsup(n, fct_solver->flux_limiter->u_tilde, fct_solver->flux_limiter->mpi_info);         printf("Paso_FCT_Solver_update_LCN: u_tilda lsup = %e (rtol = %e, max. sweeps = %d)\n",norm_u_tilde,RTOL * norm_u_tilde ,sweep_max);
        printf("Paso_FCT_Solver_update_LCN: u_tilda lsup = %e (rtol = %e, max. sweeps = %d)\n",norm_u_tilde,RTOL,sweep_max);  
262     }     }
263     errorCode = Paso_Preconditioner_Smoother_solve_byTolerance( iteration_matrix,  ((Paso_Preconditioner*) (iteration_matrix->solver_p))->gs,     errorCode = Paso_Preconditioner_Smoother_solve_byTolerance( iteration_matrix,  ((Paso_Preconditioner*) (iteration_matrix->solver_p))->gs,
264                             u, b, RTOL, &sweep_max, TRUE);                             u, b, RTOL, &sweep_max, TRUE);
# Line 311  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 319  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
319      }              }        
320      Paso_FCT_FluxLimiter_setU_tilda(flux_limiter, b); /* u_tilda = m^{-1} b */      Paso_FCT_FluxLimiter_setU_tilda(flux_limiter, b); /* u_tilda = m^{-1} b */
321      /* u_tilda_connector is completed */      /* u_tilda_connector is completed */
  for (i = 0; i < n; ++i) printf("%d : %e %e\n",i, fctp->lumped_mass_matrix[i], flux_limiter->u_tilde[i]);  
322      /**********************************************************************************************************************/        /**********************************************************************************************************************/  
323      /* calculate stopping criterium */      /* calculate stopping criterium */
324      norm_u_tilde=Paso_lsup(n, flux_limiter->u_tilde, flux_limiter->mpi_info);      norm_u_tilde=Paso_lsup(n, flux_limiter->u_tilde, flux_limiter->mpi_info);
# Line 352  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 359  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
359        
360    
361        Paso_Update(n,-1.,z,1.,b);  /* z=b-z */        Paso_Update(n,-1.,z,1.,b);  /* z=b-z */
  for (i = 0; i < n; ++i) printf("%d : %e\n",i,z[i]);  
362    
363        /* z_i += sum_{j} limitation factor_{ij} * antidiffusive_flux_{ij} */        /* z_i += sum_{j} limitation factor_{ij} * antidiffusive_flux_{ij} */
364        Paso_FCT_FluxLimiter_addLimitedFluxes_Complete(flux_limiter, z);        Paso_FCT_FluxLimiter_addLimitedFluxes_Complete(flux_limiter, z);
# Line 382  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 388  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
388                        du, z, 1, FALSE);                        du, z, 1, FALSE);
389            options->num_iter++;            options->num_iter++;
390         }                   }          
391   for (i = 0; i < n; ++i) printf("%d : %e\n",i,du[i]);          
392    
393         Paso_Update(n,1.,u,fct_solver->omega,du);         Paso_Update(n,1.,u,fct_solver->omega,du);
394         norm_du_old=norm_du;         norm_du_old=norm_du;
# Line 700  void Paso_FCT_Solver_setMuPaLu(double* o Line 706  void Paso_FCT_Solver_setMuPaLu(double* o
706    if (ABS(a)>0) {    if (ABS(a)>0) {
707        #pragma omp parallel for schedule(static) private(i, iptr_ij)        #pragma omp parallel for schedule(static) private(i, iptr_ij)
708        for (i = 0; i < n; ++i) {        for (i = 0; i < n; ++i) {
709        if ( M[i] > 0.) {
710            double sum=0;            double sum=0;
711            const double u_i=u[i];            const double u_i=u[i];
712            #pragma ivdep            #pragma ivdep
# Line 716  void Paso_FCT_Solver_setMuPaLu(double* o Line 723  void Paso_FCT_Solver_setMuPaLu(double* o
723                 sum+=l_ij*(remote_u[j]-u_i);                 sum+=l_ij*(remote_u[j]-u_i);
724            }            }
725            out[i]+=a*sum;            out[i]+=a*sum;
726        }
727        }        }
728    }    }
729  }  }

Legend:
Removed from v.3822  
changed lines
  Added in v.3824

  ViewVC Help
Powered by ViewVC 1.1.26