/[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 3795 by caltinay, Thu Feb 2 05:54:34 2012 UTC revision 3807 by gross, Mon Feb 6 06:14:25 2012 UTC
# Line 186  void Paso_FCT_Solver_initialize(const do Line 186  void Paso_FCT_Solver_initialize(const do
186  /* options2.preconditioner = PASO_GS; */  /* options2.preconditioner = PASO_GS; */
187      }      }
188      options2.use_local_preconditioner = FALSE;      options2.use_local_preconditioner = FALSE;
189      options2.sweeps=1;      options2.sweeps=-1;
190            
191      Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);      Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
192      Paso_SystemMatrix_setPreconditioner(fctp->iteration_matrix, &options2);      Paso_SystemMatrix_setPreconditioner(fctp->iteration_matrix, &options2);
# Line 226  err_t Paso_FCT_Solver_update_LCN(Paso_FC Line 226  err_t Paso_FCT_Solver_update_LCN(Paso_FC
226      Paso_SystemMatrix * iteration_matrix = fct_solver->transportproblem->iteration_matrix;      Paso_SystemMatrix * iteration_matrix = fct_solver->transportproblem->iteration_matrix;
227      err_t errorCode = SOLVER_NO_ERROR;      err_t errorCode = SOLVER_NO_ERROR;
228        
     
229      Paso_Coupler_startCollect(fct_solver->u_old_coupler,u_old);      Paso_Coupler_startCollect(fct_solver->u_old_coupler,u_old);
230      Paso_Coupler_finishCollect(fct_solver->u_old_coupler);      Paso_Coupler_finishCollect(fct_solver->u_old_coupler);
231            
# Line 235  err_t Paso_FCT_Solver_update_LCN(Paso_FC Line 234  err_t Paso_FCT_Solver_update_LCN(Paso_FC
234          * low order transport matrix l. Therefore a=-dt*0.5 is used. */          * low order transport matrix l. Therefore a=-dt*0.5 is used. */
235                
236      Paso_FCT_Solver_setMuPaLu(b, fct_solver->transportproblem->lumped_mass_matrix,      Paso_FCT_Solver_setMuPaLu(b, fct_solver->transportproblem->lumped_mass_matrix,
237                   fct_solver->u_old_coupler, -dt*0.5, iteration_matrix);                   fct_solver->u_old_coupler, -dt*0.5, iteration_matrix);
       
238      /* solve for u_tilde : u_tilda = m^{-1} * b   */      /* solve for u_tilde : u_tilda = m^{-1} * b   */
239      Paso_FCT_FluxLimiter_setU_tilda(fct_solver->flux_limiter, b);      Paso_FCT_FluxLimiter_setU_tilda(fct_solver->flux_limiter, b);
240      /* u_tilda_connector is completed */      /* u_tilda_connector is completed */
241          
242      /* calculate anti-diffusive fluxes for u_tilda */      /* calculate anti-diffusive fluxes for u_tilda */
243      Paso_FCT_setAntiDiffusionFlux_linearCN(fct_solver->flux_limiter->antidiffusive_fluxes,      Paso_FCT_setAntiDiffusionFlux_linearCN(fct_solver->flux_limiter->antidiffusive_fluxes,
244                         fct_solver->transportproblem, dt,                         fct_solver->transportproblem, dt,
245                         fct_solver->flux_limiter->u_tilde_coupler,                         fct_solver->flux_limiter->u_tilde_coupler,
246                         fct_solver->u_old_coupler);                         fct_solver->u_old_coupler);
247      
248    
249     /* b_i += sum_{j} limitation factor_{ij} * antidiffusive_flux_{ij} */     /* b_i += sum_{j} limitation factor_{ij} * antidiffusive_flux_{ij} */
250     Paso_FCT_FluxLimiter_addLimitedFluxes_Start(fct_solver->flux_limiter);     Paso_FCT_FluxLimiter_addLimitedFluxes_Start(fct_solver->flux_limiter);
251     Paso_FCT_FluxLimiter_addLimitedFluxes_Complete(fct_solver->flux_limiter, b);     Paso_FCT_FluxLimiter_addLimitedFluxes_Complete(fct_solver->flux_limiter, b);
252        
253        
254     /* solve (m-dt/2*L) u = b */       /* solve (m-dt/2*L) u = b in the form (omega*m-L) u = b * omega with omega*dt/2=1 */
255      
256     /* initial guess is u<- -u + 2*u_tilde */     /* initial guess is u<- -u + 2*u_tilde */
257     Paso_Update(n, -1., u, 2., fct_solver->flux_limiter->u_tilde);     Paso_Update(n, -1., u, 2., fct_solver->flux_limiter->u_tilde);
258       Paso_Scale(n, b,fct_solver->omega );
259     sweep_max = MAX((int) (- 2 * log(RTOL)/log(2.)-0.5),1);     sweep_max = MAX((int) (- 2 * log(RTOL)/log(2.)-0.5),1);
260      
261       if (options->verbose) {
262           const double norm_u_tilde=Paso_lsup(n, fct_solver->flux_limiter->u_tilde, fct_solver->flux_limiter->mpi_info);
263           printf("Paso_FCT_Solver_update_LCN: u_tilda lsup = %e (rtol = %e, max. sweeps = %d)\n",norm_u_tilde,RTOL,sweep_max);
264       }
265     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,
266                             u, b, RTOL, &sweep_max, TRUE);                             u, b, RTOL, &sweep_max, TRUE);
267     if (errorCode == PRECONDITIONER_NO_ERROR) {     if (errorCode == PRECONDITIONER_NO_ERROR) {
# Line 551  void Paso_FCT_setAntiDiffusionFlux_linea Line 556  void Paso_FCT_setAntiDiffusionFlux_linea
556    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
557        const double u_tilde_i = u_tilde[i];        const double u_tilde_i = u_tilde[i];
558        const double u_old_i   = u_old[i];        const double u_old_i   = u_old[i];
559        const double du_i      = u_old_i - u_tilde_i;        const double du_i      = u_tilde_i - u_old_i;
560        #pragma ivdep        #pragma ivdep
561        for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {        for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
562            
# Line 560  void Paso_FCT_setAntiDiffusionFlux_linea Line 565  void Paso_FCT_setAntiDiffusionFlux_linea
565        const double d_ij    = fct->transport_matrix->mainBlock->val[iptr_ij]+fct->iteration_matrix->mainBlock->val[iptr_ij]; /* this is in fact -d_ij */        const double d_ij    = fct->transport_matrix->mainBlock->val[iptr_ij]+fct->iteration_matrix->mainBlock->val[iptr_ij]; /* this is in fact -d_ij */
566            const double u_tilde_j = u_tilde[j];            const double u_tilde_j = u_tilde[j];
567        const double u_old_j = u_old[j];        const double u_old_j = u_old[j];
568        const double du_j    = u_old_j - u_tilde_j;        const double du_j    = u_tilde_j - u_old_j;
569                
570        flux_matrix->mainBlock->val[iptr_ij]=2 * m_ij * ( du_i - du_j ) + dt * d_ij * ( u_tilde_j - u_tilde_i);        flux_matrix->mainBlock->val[iptr_ij]= 2 * m_ij * ( du_i - du_j ) - dt * d_ij * ( u_tilde_i - u_tilde_j);
   
571        }        }
572        #pragma ivdep        #pragma ivdep
573        for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {        for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
# Line 573  void Paso_FCT_setAntiDiffusionFlux_linea Line 577  void Paso_FCT_setAntiDiffusionFlux_linea
577      const double d_ij    = fct->transport_matrix->col_coupleBlock->val[iptr_ij]+fct->iteration_matrix->col_coupleBlock->val[iptr_ij];/* this is in fact -d_ij */      const double d_ij    = fct->transport_matrix->col_coupleBlock->val[iptr_ij]+fct->iteration_matrix->col_coupleBlock->val[iptr_ij];/* this is in fact -d_ij */
578          const double u_tilde_j = remote_u_tilde[j];          const double u_tilde_j = remote_u_tilde[j];
579      const double u_old_j = remote_u_old[j];      const double u_old_j = remote_u_old[j];
580      const double du_j    = u_old_j - u_tilde_j;      const double du_j    = u_tilde_j - u_old_j;
581    
582      flux_matrix->col_coupleBlock->val[iptr_ij]= 2 * m_ij * ( du_i - du_j ) + dt * d_ij *  ( u_tilde_j - u_tilde_i);      flux_matrix->col_coupleBlock->val[iptr_ij]= 2 * m_ij * ( du_i - du_j ) - dt * d_ij *  ( u_tilde_i - u_tilde_j);
583            
584        }        }
585    }    }

Legend:
Removed from v.3795  
changed lines
  Added in v.3807

  ViewVC Help
Powered by ViewVC 1.1.26