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

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

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

revision 1371 by gross, Thu Jan 3 06:11:21 2008 UTC revision 1372 by gross, Thu Jan 3 06:30:47 2008 UTC
# Line 89  void Paso_SolverFCT_solve(Paso_FCTranspo Line 89  void Paso_SolverFCT_solve(Paso_FCTranspo
89     index_t i;     index_t i;
90     int n_substeps,n;     int n_substeps,n;
91     double dt2=fctp->dt_max, dt2_loc, rtmp,rtmp2,t;     double dt2=fctp->dt_max, dt2_loc, rtmp,rtmp2,t;
92       dim_t n_rows=Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);
93     if (dt<=0.) {     if (dt<=0.) {
94         Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");         Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");
95     }     }
# Line 100  void Paso_SolverFCT_solve(Paso_FCTranspo Line 101  void Paso_SolverFCT_solve(Paso_FCTranspo
101          Paso_FCTransportProblem_addAdvectivePart(fctp,1.);          Paso_FCTransportProblem_addAdvectivePart(fctp,1.);
102          /* create a copy of the main diagonal entires of the transport-matrix */          /* create a copy of the main diagonal entires of the transport-matrix */
103          #pragma omp parallel for schedule(static) private(i)          #pragma omp parallel for schedule(static) private(i)
104          for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {          for (i=0;i<n_rows;++i) {
105                 fctp->transport_matrix_diagonal[i]=                 fctp->transport_matrix_diagonal[i]=
106                          fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];                          fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];
107          }          }
# Line 118  Paso_SystemMatrix_saveMM(fctp->transport Line 119  Paso_SystemMatrix_saveMM(fctp->transport
119              {              {
120                 dt2_loc=LARGE_POSITIVE_FLOAT;                 dt2_loc=LARGE_POSITIVE_FLOAT;
121                 #pragma omp for schedule(static) private(i,rtmp,rtmp2)                 #pragma omp for schedule(static) private(i,rtmp,rtmp2)
122                 for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {                 for (i=0;i<n_rows;++i) {
123                      rtmp=fctp->transport_matrix_diagonal[i];                      rtmp=fctp->transport_matrix_diagonal[i];
124                      rtmp2=fctp->lumped_mass_matrix[i];                      rtmp2=fctp->lumped_mass_matrix[i];
125                      if ( (rtmp<0 && rtmp2>=0.) || (rtmp>0 && rtmp2<=0.) ) {                      if ( (rtmp<0 && rtmp2>=0.) || (rtmp>0 && rtmp2<=0.) ) {
126                          dt2_loc=MIN(dt2_loc,-rtmp2/rtmp);                          dt2_loc=MIN(dt2_loc,-rtmp2/rtmp);
127                      }                      }
128                  }                  }
129                  #pragma omp critcal                  #pragma omp critical
130                  {                  {
131                      dt2=MIN(dt2,dt2_loc);                      dt2=MIN(dt2,dt2_loc);
132                  }                  }
# Line 169  Paso_SystemMatrix_saveMM(fctp->transport Line 170  Paso_SystemMatrix_saveMM(fctp->transport
170    
171     if (fctp->theta>0) {     if (fctp->theta>0) {
172          #pragma omp parallel for schedule(static) private(i)          #pragma omp parallel for schedule(static) private(i)
173          for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {          for (i=0;i<n_rows;++i) {
174                 fctp->transport_matrix_diagonal[i]=                 fctp->transport_matrix_diagonal[i]=
175                          fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];                          fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];
176          }          }
# Line 182  Paso_SystemMatrix_saveMM(fctp->transport Line 183  Paso_SystemMatrix_saveMM(fctp->transport
183          printf("substep step %d at t=%e\n",n+1,t);          printf("substep step %d at t=%e\n",n+1,t);
184          Paso_FCTransportProblem_setFlux(fctp,fctp->u,u); /* u stores F(u_last)*/          Paso_FCTransportProblem_setFlux(fctp,fctp->u,u); /* u stores F(u_last)*/
185          #pragma omp parallel for schedule(static) private(i)          #pragma omp parallel for schedule(static) private(i)
186          for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {          for (i=0;i<n_rows;++i) {
187              rtmp=fctp->u[i];              rtmp=fctp->u[i];
188              fctp->u[i]=fctp->u[i]+dt2*u[i]/fctp->lumped_mass_matrix[i];              fctp->u[i]=fctp->u[i]+dt2*u[i]/fctp->lumped_mass_matrix[i];
189              printf("%d : %e %e %e\n",i,u[i],fctp->u[i],rtmp);              printf("%d : %e %e %e\n",i,u[i],fctp->u[i],rtmp);
# Line 193  Paso_SystemMatrix_saveMM(fctp->transport Line 194  Paso_SystemMatrix_saveMM(fctp->transport
194    }    }
195    if (Paso_noError()) {    if (Paso_noError()) {
196       #pragma omp parallel for schedule(static) private(i)       #pragma omp parallel for schedule(static) private(i)
197       for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {       for (i=0;i<n_rows;++i) {
198         u[i]=fctp->u[i];         u[i]=fctp->u[i];
199       }       }
200    }    }
Legend:
Removed from v.1371  
changed lines
  Added in v.1372
  ViewVC Help
Powered by ViewVC 1.1.26