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

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

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

revision 3004 by gross, Tue Mar 16 01:32:43 2010 UTC revision 3005 by gross, Thu Apr 22 05:59:31 2010 UTC
# Line 54  void Paso_ReactiveSolver_free(Paso_React Line 54  void Paso_ReactiveSolver_free(Paso_React
54  }  }
55  double Paso_ReactiveSolver_getSafeTimeStepSize(Paso_TransportProblem* fctp)  double Paso_ReactiveSolver_getSafeTimeStepSize(Paso_TransportProblem* fctp)
56  {  {
57     return LARGE_POSITIVE_FLOAT;    double dt_max=LARGE_POSITIVE_FLOAT;
58    
59    
60    
61      double beta_0=0., beta_0_loc, betap, beta=EPSILON;
62      dim_t i, n;
63      int fail=0;
64      double dt_max_loc;
65      register double m_ii,m_i, d_ii;
66      #ifdef PASO_MPI
67         double rtmp[2], rtmp_loc[2];
68      #endif
69      n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
70    
71      #pragma omp parallel private(beta_0_loc)
72      {
73          beta_0_loc=0.;
74          #pragma omp for schedule(static) private(i,m_ii,m_i, d_ii)
75          for (i=0;i<n;++i) {
76             m_i=fctp->lumped_mass_matrix[i];
77             m_ii=fctp->main_diagonal_mass_matrix[i];
78    /* printf(" %d : %e %e -> %e\n", i, m_i, m_ii, m_i/m_ii); */
79             if ( m_ii>0 ) {
80                beta_0_loc=MAX(m_i/m_ii, beta_0_loc);
81             } else {
82                fail=1;
83             }
84          }
85          #pragma omp critical
86          {
87              beta_0=MAX(beta_0,beta_0_loc);
88          }
89      }
90      #ifdef PASO_MPI
91         rtmp_loc[0]=beta_0;
92         rtmp_loc[1]=1.(*(double) fail );
93         MPI_Allreduce(rtmp_loc, rtmp, 2, MPI_DOUBLE, MPI_MAX, fctp->mpi_info->comm);
94         beta_0=rtmp_loc[0];
95         fail= (rtmp_loc[0] > 0) ? 1 : 0;
96      #endif
97      if (fail>0) {
98         Paso_setError(TYPE_ERROR,"Paso_ReactiveSolver_getSafeTimeStepSize: non-positive main diagonal entry in mass matrix.");
99      }
100      
101      if ( Paso_noError()) {
102          if ( ( beta_0 < 2. ) && ( beta_0 >= 1. ) ) {
103              beta=MAX(2*(beta_0-1.)/beta_0, EPSILON);
104           } else {
105              Paso_setError(TYPE_ERROR,"Paso_ReactiveSolver_getSafeTimeStepSize: non-positive main diagonal entry in mass matrix.");
106           }
107      }
108    printf("beta = %e from beta_0 = %e\n",beta, beta_0-1);
109      if ( Paso_noError()) {
110        
111         #pragma omp parallel private(dt_max_loc, betap)
112         {
113             betap=beta+1.;
114             dt_max_loc=LARGE_POSITIVE_FLOAT;
115             #pragma omp for schedule(static) private(i, m_ii, m_i, d_ii)
116             for (i=0;i<n;++i) {
117                 d_ii=fctp->reactive_matrix[i];
118                 m_i=fctp->lumped_mass_matrix[i];
119                 m_ii=fctp->main_diagonal_mass_matrix[i];
120    
121                 if (d_ii > 0.)  dt_max_loc=MIN(dt_max_loc, (betap* m_ii - m_i)/d_ii );
122             }
123             #pragma omp critical
124             {
125                 dt_max=MIN(dt_max,dt_max_loc);
126            }
127         }
128         dt_max=1./(beta* Paso_Transport_getTheta(fctp));
129          #ifdef PASO_MPI
130             dt_max_loc=dt_max;
131             MPI_Allreduce(&dt_max_loc, &dt_max, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
132          #endif
133      }
134      return dt_max;
135  }  }

Legend:
Removed from v.3004  
changed lines
  Added in v.3005

  ViewVC Help
Powered by ViewVC 1.1.26