/[escript]/branches/symbolic_from_3470/paso/src/ReactiveSolver.c
ViewVC logotype

Diff of /branches/symbolic_from_3470/paso/src/ReactiveSolver.c

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

revision 3867 by caltinay, Thu Feb 9 00:27:46 2012 UTC revision 3868 by caltinay, Thu Mar 15 06:07:08 2012 UTC
# Line 41  err_t  Paso_ReactiveSolver_solve(Paso_Re Line 41  err_t  Paso_ReactiveSolver_solve(Paso_Re
41       const double EXP_LIM_MAX =PASO_RT_EXP_LIM_MAX;       const double EXP_LIM_MAX =PASO_RT_EXP_LIM_MAX;
42       const double dt = support->dt;       const double dt = support->dt;
43       index_t fail=0;       index_t fail=0;
44       register double d_ii, m_i, x_i, e_i, u_i, F_i;       register double d_ii, x_i, e_i, u_i, F_i;
45       dim_t i;       dim_t i;
46       const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);       const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
47            
48       #pragma omp parallel for schedule(static) private(i, d_ii, m_i, x_i, e_i, u_i, F_i)       #pragma omp parallel for schedule(static) private(i, d_ii, x_i, e_i, u_i, F_i)
49       for (i=0;i<n;++i) {       for (i=0;i<n;++i) {
50          d_ii=fctp->reactive_matrix[i];          const double  m_i=fctp->lumped_mass_matrix[i];
51          m_i=fctp->lumped_mass_matrix[i];          if (m_i>0) {
52          x_i=dt*d_ii/m_i;          
53      if (x_i >= EXP_LIM_MAX) {          d_ii=fctp->reactive_matrix[i];
54        fail=1;          x_i=dt*d_ii/m_i;
55      } else  {          if (x_i >= EXP_LIM_MAX) {
56          F_i=source[i];            fail=1;
57          e_i=exp(x_i);          } else  {
58          u_i=e_i*u_old[i];          F_i=source[i];
59          if ( abs(x_i) > EXP_LIM_MIN) {          e_i=exp(x_i);
60          u_i+=F_i/d_ii*(e_i-1.);          u_i=e_i*u_old[i];
61          } else {          if ( abs(x_i) > EXP_LIM_MIN) {
62          u_i+=F_i*dt/m_i * (1. + x_i/2); /* second order approximation of ( exp(x_i)-1)/x_i */              u_i+=F_i/d_ii*(e_i-1.);
63            } else {
64                u_i+=F_i*dt/m_i * (1. + x_i/2); /* second order approximation of ( exp(x_i)-1)/x_i */
65            }
66            u[i]=u_i;
67          }          }
68          u[i]=u_i;      } else {
69            u[i]=u_old[i] + dt * source[i] ; /* constraints added */
70      }      }
71      }      }
72      #ifdef ESYS_MPI      #ifdef ESYS_MPI
# Line 95  double Paso_ReactiveSolver_getSafeTimeSt Line 100  double Paso_ReactiveSolver_getSafeTimeSt
100  {  {
101    
102       const double EXP_LIM_MAX =PASO_RT_EXP_LIM_MAX;       const double EXP_LIM_MAX =PASO_RT_EXP_LIM_MAX;
103        
104       double dt_max=LARGE_POSITIVE_FLOAT, dt_max_loc;         double dt_max=LARGE_POSITIVE_FLOAT, dt_max_loc;  
      index_t fail_loc, fail;  
105       dim_t i;       dim_t i;
106       const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);       const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
107       register double d_ii,m_i;       register double d_ii,m_i;
# Line 105  double Paso_ReactiveSolver_getSafeTimeSt Line 110  double Paso_ReactiveSolver_getSafeTimeSt
110           *  calculate time step size:                                                     *  calculate time step size:                                          
111          */          */
112          dt_max=LARGE_POSITIVE_FLOAT;          dt_max=LARGE_POSITIVE_FLOAT;
113      fail=0;          #pragma omp parallel private(dt_max_loc)
         #pragma omp parallel private(dt_max_loc, fail_loc)  
114          {          {
115                 dt_max_loc=LARGE_POSITIVE_FLOAT;                  dt_max_loc=LARGE_POSITIVE_FLOAT;
            fail_loc=0;  
116                 #pragma omp for schedule(static) private(i,d_ii,m_i)                 #pragma omp for schedule(static) private(i,d_ii,m_i)
117                 for (i=0;i<n;++i) {                 for (i=0;i<n;++i) {
118                    d_ii=fctp->reactive_matrix[i];                    d_ii=fctp->reactive_matrix[i];
119                    m_i=fctp->lumped_mass_matrix[i];                    m_i=fctp->lumped_mass_matrix[i];
120            if (m_i > 0) {            if (m_i > 0) { /* no constraint */
121                if ( d_ii>0 ) dt_max_loc=MIN(dt_max_loc, m_i/d_ii);                if ( d_ii>0 ) dt_max_loc=MIN(dt_max_loc, m_i/d_ii);
           } else {  
               fail_loc=-1;  
122            }            }
123                 }                 }
124                 #pragma omp critical                 #pragma omp critical
125                 {                 {
126                    dt_max=MIN(dt_max, dt_max_loc);                    dt_max=MIN(dt_max, dt_max_loc);
           fail=MIN(fail, fail_loc);  
127                 }                 }
128          }          }
129          #ifdef ESYS_MPI          #ifdef ESYS_MPI
130          {          {
131             double rtmp_loc[2], rtmp[2];             dt_max_loc=dt_max;
132                 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;  
133      }      }
134          #endif          #endif
135          if (fail < 0 ) {  
136         Esys_setError(VALUE_ERROR, "Paso_ReactiveSolver_getSafeTimeStepSize: negative mass term detected.");      if (dt_max < LARGE_POSITIVE_FLOAT ) {
137         return -1;              dt_max*=0.5*EXP_LIM_MAX;                /* make sure there is no exp overflow */
138      } else {          } else {
         if (dt_max < LARGE_POSITIVE_FLOAT ) {  
                dt_max*=0.5*EXP_LIM_MAX;  
             } else {  
139                 dt_max=LARGE_POSITIVE_FLOAT;                 dt_max=LARGE_POSITIVE_FLOAT;
140              }          }
     }  
141     }     }
142     return dt_max;     return dt_max;
143  }  }

Legend:
Removed from v.3867  
changed lines
  Added in v.3868

  ViewVC Help
Powered by ViewVC 1.1.26