/[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 3807 by gross, Mon Feb 6 06:14:25 2012 UTC revision 3822 by gross, Mon Feb 13 06:22:06 2012 UTC
# Line 98  double Paso_FCT_Solver_getSafeTimeStepSi Line 98  double Paso_FCT_Solver_getSafeTimeStepSi
98  {  {
99     dim_t i, n;     dim_t i, n;
100     double dt_max=LARGE_POSITIVE_FLOAT;     double dt_max=LARGE_POSITIVE_FLOAT;
    index_t fail=0;  
101     n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
102     /* set low order transport operator */     /* set low order transport operator */
103     Paso_FCT_setLowOrderOperator(fctp);     Paso_FCT_setLowOrderOperator(fctp);
# Line 108  double Paso_FCT_Solver_getSafeTimeStepSi Line 107  double Paso_FCT_Solver_getSafeTimeStepSi
107           *  calculate time step size:                                                     *  calculate time step size:                                          
108          */          */
109          dt_max=LARGE_POSITIVE_FLOAT;          dt_max=LARGE_POSITIVE_FLOAT;
     fail=0;  
110          #pragma omp parallel private(i)          #pragma omp parallel private(i)
111          {          {
112                 double dt_max_loc=LARGE_POSITIVE_FLOAT;                 double dt_max_loc=LARGE_POSITIVE_FLOAT;
            index_t fail_loc=0;  
113                 #pragma omp for schedule(static)                 #pragma omp for schedule(static)
114                 for (i=0;i<n;++i) {                 for (i=0;i<n;++i) {
115                    const double l_ii=fctp->main_diagonal_low_order_transport_matrix[i];                    const double l_ii=fctp->main_diagonal_low_order_transport_matrix[i];
116                    const double m_i=fctp->lumped_mass_matrix[i];                    const double m_i=fctp->lumped_mass_matrix[i];
117            if ( (m_i > 0) ) {            if ( m_i > 0 ) {
118                if (l_ii<0) dt_max_loc=MIN(dt_max_loc,m_i/(-l_ii));                if (l_ii<0) dt_max_loc=MIN(dt_max_loc,m_i/(-l_ii));
119            } else {                    }
120                fail_loc=-1;             }
           }  
                }  
121                 #pragma omp critical                 #pragma omp critical
122                 {                 {
123                    dt_max=MIN(dt_max,dt_max_loc);                    dt_max=MIN(dt_max,dt_max_loc);
           fail=MIN(fail, fail_loc);  
124                 }                 }
125          }          }
126          #ifdef ESYS_MPI          #ifdef ESYS_MPI
127          {          {
128             double rtmp_loc[2], rtmp[2];                 dt_max_loc=dt_max;
129                 rtmp_loc[0]=dt_max;                 MPI_Allreduce(&dt_max_loc, &dt_max, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
            rtmp_loc[1]= (double) fail;  
                MPI_Allreduce(rtmp_loc, rtmp, 2, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);  
            dt_max=rtmp[0];  
            fail = rtmp[1] < 0 ? -1 : 0;  
130      }      }
131          #endif          #endif
132          if (fail < 0 ) {      if (dt_max<LARGE_POSITIVE_FLOAT) dt_max*=2.;
        Esys_setError(VALUE_ERROR, "Paso_FCTSolver_getSafeTimeStepSize: negative mass matrix entries detected.");  
        return -1;  
     } else {  
         if (dt_max<LARGE_POSITIVE_FLOAT) dt_max*=2.;  
     }  
133     }     }
134     return dt_max;     return dt_max;
135  }  }
# Line 171  void Paso_FCT_Solver_initialize(const do Line 156  void Paso_FCT_Solver_initialize(const do
156      fct_solver->dt = dt;      fct_solver->dt = dt;
157      #pragma omp parallel for private(i)      #pragma omp parallel for private(i)
158      for (i = 0; i < n; ++i) {      for (i = 0; i < n; ++i) {
159             const double m=fctp->lumped_mass_matrix[i];             const double m_i=fctp->lumped_mass_matrix[i];
160         const double l_ii = fctp->main_diagonal_low_order_transport_matrix[i];         const double l_ii = fctp->main_diagonal_low_order_transport_matrix[i];
161             fctp->iteration_matrix->mainBlock->val[main_iptr[i]] = m * omega - l_ii;         if ( m_i > 0 ) {
162                   fctp->iteration_matrix->mainBlock->val[main_iptr[i]] = m_i * omega - l_ii;
163           } else {
164               fctp->iteration_matrix->mainBlock->val[main_iptr[i]] = ABS(m_i * omega - l_ii)/(EPSILON*EPSILON);
165           }
166      }      }
167            
168      /* allocate preconditioner/solver */      /* allocate preconditioner/solver */
# Line 229  err_t Paso_FCT_Solver_update_LCN(Paso_FC Line 218  err_t Paso_FCT_Solver_update_LCN(Paso_FC
218      Paso_Coupler_startCollect(fct_solver->u_old_coupler,u_old);      Paso_Coupler_startCollect(fct_solver->u_old_coupler,u_old);
219      Paso_Coupler_finishCollect(fct_solver->u_old_coupler);      Paso_Coupler_finishCollect(fct_solver->u_old_coupler);
220            
221      /* b[i]=m*u_tilde[i] = m u_old[i] + dt/2 sum_{j <> i} l_{ij}*(u_old[j]-u_old[i])      /* b[i]=m*u_tilde[i] = m u_old[i] + dt/2 sum_{j <> i} l_{ij}*(u_old[j]-u_old[i])  
222               = u_tilde[i]   = u_old[i] where constraint m<0.
223          * note that iteration_matrix stores the negative values of the          * note that iteration_matrix stores the negative values of the
224          * low order transport matrix l. Therefore a=-dt*0.5 is used. */          * low order transport matrix l. Therefore a=-dt*0.5 is used. */
225                
# Line 306  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 296  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
296            /* b[i]=m_i* u_old[i] */            /* b[i]=m_i* u_old[i] */
297            #pragma omp for private(i) schedule(static)            #pragma omp for private(i) schedule(static)
298            for (i = 0; i < n; ++i) {            for (i = 0; i < n; ++i) {
299                 b[i]=u_old[i]* fctp->lumped_mass_matrix[i];             if (fctp->lumped_mass_matrix[i] > 0 ) {
300               b[i]=u_old[i]* fctp->lumped_mass_matrix[i];
301               } else {
302               b[i]=u_old[i];
303               }
304            }            }
305      } else {      } else {
306         /* b[i]=m_i* u_old[i] + dt/2 sum_{j <> i} l_{ij}*(u_old[j]-u_old[i]) = m_i * u_tilde_i         /* b[i]=m_i* u_old[i] + dt/2 sum_{j <> i} l_{ij}*(u_old[j]-u_old[i]) = m_i * u_tilde_i where m_i>0
307        *     = u_old[i]  otherwise
308          * note that iteration_matrix stores the negative values of the          * note that iteration_matrix stores the negative values of the
309          * low order transport matrix l. Therefore a=-dt*0.5 is used. */          * low order transport matrix l. Therefore a=-dt*0.5 is used. */
310         Paso_FCT_Solver_setMuPaLu(b,fctp->lumped_mass_matrix,fct_solver->u_old_coupler,-dt*0.5,fctp->iteration_matrix);         Paso_FCT_Solver_setMuPaLu(b,fctp->lumped_mass_matrix,fct_solver->u_old_coupler,-dt*0.5,fctp->iteration_matrix);
311      }              }        
312      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 */
313      /* u_tilda_connector is completed */      /* u_tilda_connector is completed */
314       for (i = 0; i < n; ++i) printf("%d : %e %e\n",i, fctp->lumped_mass_matrix[i], flux_limiter->u_tilde[i]);
315      /**********************************************************************************************************************/        /**********************************************************************************************************************/  
316      /* calculate stopping criterium */      /* calculate stopping criterium */
317      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 341  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 336  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
336           Paso_FCT_FluxLimiter_addLimitedFluxes_Start(flux_limiter); /* uses u_tilde */             Paso_FCT_FluxLimiter_addLimitedFluxes_Start(flux_limiter); /* uses u_tilde */  
337    
338           /*           /*
339        * z_m[i]=b[i] - (m_i*u[i] - omega*sum_{j<>i} l_{ij} (u[j]-u[i]) ) omega = dt/2 or dt .        * z_m[i]=b[i] - (m_i*u[i] - omega*sum_{j<>i} l_{ij} (u[j]-u[i]) ) where m_i>0
340          *       ==b[i] - u[i] = u_tilda[i]-u[i] =0 otherwise
341          *
342          * omega = dt/2 or dt .
343        *        *
344        * note that iteration_matrix stores the negative values of the        * note that iteration_matrix stores the negative values of the
345        * low order transport matrix l. Therefore a=dt*theta is used.        * low order transport matrix l. Therefore a=dt*theta is used.
# Line 354  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 352  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
352        
353    
354        Paso_Update(n,-1.,z,1.,b);  /* z=b-z */        Paso_Update(n,-1.,z,1.,b);  /* z=b-z */
355     for (i = 0; i < n; ++i) printf("%d : %e\n",i,z[i]);
356    
357        /* z_i += sum_{j} limitation factor_{ij} * antidiffusive_flux_{ij} */        /* z_i += sum_{j} limitation factor_{ij} * antidiffusive_flux_{ij} */
358        Paso_FCT_FluxLimiter_addLimitedFluxes_Complete(flux_limiter, z);        Paso_FCT_FluxLimiter_addLimitedFluxes_Complete(flux_limiter, z);
# Line 384  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 382  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
382                        du, z, 1, FALSE);                        du, z, 1, FALSE);
383            options->num_iter++;            options->num_iter++;
384         }                   }          
385           for (i = 0; i < n; ++i) printf("%d : %e\n",i,du[i]);      
386    
387         Paso_Update(n,1.,u,fct_solver->omega,du);         Paso_Update(n,1.,u,fct_solver->omega,du);
388         norm_du_old=norm_du;         norm_du_old=norm_du;
# Line 674  printf("a[%d,%d]=%e\n",j,i,rtmp2); Line 672  printf("a[%d,%d]=%e\n",j,i,rtmp2);
672  }  }
673    
674  /*  /*
675   * out_i=m_i u_i + a * \sum_{j <> i} l_{ij} (u_j-u_i)   * out_i=m_i u_i + a * \sum_{j <> i} l_{ij} (u_j-u_i) where m_i>0
676     *       = u_i                                        where m_i<=0
677   *   *
678   */   */
679  void Paso_FCT_Solver_setMuPaLu(double* out,  void Paso_FCT_Solver_setMuPaLu(double* out,
# Line 692  void Paso_FCT_Solver_setMuPaLu(double* o Line 691  void Paso_FCT_Solver_setMuPaLu(double* o
691    
692    #pragma omp parallel for private(i) schedule(static)    #pragma omp parallel for private(i) schedule(static)
693    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
694        out[i]=M[i]*u[i];        if ( M[i] > 0.) {
695              out[i]=M[i]*u[i];
696          } else {
697          out[i]=u[i];
698          }
699    }    }
700    if (ABS(a)>0) {    if (ABS(a)>0) {
701        #pragma omp parallel for schedule(static) private(i, iptr_ij)        #pragma omp parallel for schedule(static) private(i, iptr_ij)

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

  ViewVC Help
Powered by ViewVC 1.1.26