/[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 1407 by gross, Mon Feb 4 06:45:48 2008 UTC revision 1410 by gross, Thu Feb 7 04:24:00 2008 UTC
# Line 154  void Paso_SolverFCT_solve(Paso_FCTranspo Line 154  void Paso_SolverFCT_solve(Paso_FCTranspo
154         }         }
155         dt2=dt/n_substeps;         dt2=dt/n_substeps;
156         printf("%d time steps of size is %e (theta = %e, dt_max=%e).\n",n_substeps, dt2,fctp->theta, dt_max);         printf("%d time steps of size is %e (theta = %e, dt_max=%e).\n",n_substeps, dt2,fctp->theta, dt_max);
   
157         /*         /*
158      * seperate source into positive and negative part:      * seperate source into positive and negative part:
159      */      */
 printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF\n");  
160          #pragma omp parallel for private(i,rtmp)          #pragma omp parallel for private(i,rtmp)
161          for (i = 0; i < n_rows; ++i) {          for (i = 0; i < n_rows; ++i) {
162            rtmp=source[i];            rtmp=source[i];
163            if (rtmp <0) {            if (rtmp <0) {
164               sourceP[i]=-rtmp;               sourceN[i]=-rtmp;
165                 sourceP[i]=0;
166            } else {            } else {
167               sourceN[i]= rtmp;               sourceN[i]= 0;
168                 sourceP[i]= rtmp;
169            }            }
170          }          }
171  printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF\n");  /*        for (i = 0; i < n_rows; ++i) printf("%d : %e \n",i,source[i],sourceP[i],sourceN[i]) */
172          /*          /*
173           * now the show can begin:           * let the show begin!!!!
174           *           *
175           */           */
         #pragma omp parallel for schedule(static) private(i)  
         for (i=0;i<n_rows;++i) u[i]=fctp->u[i]-fctp->u_min;  
 printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF\n");  
   
176          t=dt2;          t=dt2;
177          n=0;          n=0;
178          tolerance=options->tolerance;          tolerance=options->tolerance;
179          while(n<n_substeps && Paso_noError()) {          while(n<n_substeps && Paso_noError()) {
 printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF\n");  
180              printf("substep step %d at t=%e\n",n+1,t);              printf("substep step %d at t=%e\n",n+1,t);
181                #pragma omp parallel for private(i)
182                 for (i = 0; i < n_rows; ++i) {
183                          u[i]=fctp->u[i];
184                 }
185              /*              /*
186               * b^n[i]=m u^n[i] + dt*(1-theta) sum_{j <> i} l_{ij}*(u^n[j]-u^n[i]) + dt*sourceP[i]               * b^n[i]=m u^n[i] + dt2*(1-theta) sum_{j <> i} l_{ij}*(u^n[j]-u^n[i]) + dt2*sourceP[i]
187               *               *
188               * note that iteration_matrix stores the negative values of the               * note that iteration_matrix stores the negative values of the
189               * low order transport matrix l therefore a=-dt*(1-fctp->theta) is used.               * low order transport matrix l therefore a=-dt2*(1-fctp->theta) is used.
190               *               *
191               */               */
192               Paso_SolverFCT_setMuPaLuPbQ(b_n,fctp->lumped_mass_matrix,u,               Paso_SolverFCT_setMuPaLuPbQ(b_n,fctp->lumped_mass_matrix,u,
193                                            -dt*(1-fctp->theta),fctp->iteration_matrix,dt,sourceP);                                            -dt2*(1-fctp->theta),fctp->iteration_matrix,dt2,sourceP);
194               /*               /*
195            *   uTilde_n[i]=b[i]/m[i]            *   uTilde_n[i]=b[i]/m[i]
196            *            *
197            *   a[i,i]=m[i]/(dt theta) + \frac{1}{\theta} \frac{q^-[i]}-l[i,i]            *   a[i,i]=m[i]/(dt2 theta) + \frac{1}{\theta} \frac{q^-[i]}-l[i,i]
198            *            *
199            */            */
 printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF\n");  
200               if (fctp->theta > 0) {               if (fctp->theta > 0) {
201                   Paso_solve_free(fctp->iteration_matrix);                   Paso_solve_free(fctp->iteration_matrix);
202                   rtmp1=1./(dt*fctp->theta);                   omega=1./(dt2*fctp->theta);
203                   rtmp2=1./fctp->theta;                   rtmp2=dt2*omega;
204                   #pragma omp parallel for private(i,rtmp,rtmp3,rtmp4)                   #pragma omp parallel for private(i,rtmp,rtmp3,rtmp4)
205                   for (i = 0; i < n_rows; ++i) {                   for (i = 0; i < n_rows; ++i) {
206                        rtmp=fctp->lumped_mass_matrix[i];                        rtmp=fctp->lumped_mass_matrix[i];
# Line 211  printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF Line 209  printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
209                        } else {                        } else {
210                           rtmp3=u[i];                           rtmp3=u[i];
211                        }                        }
212                        rtmp4=rtmp*rtmp1-fctp->main_diagonal_low_order_transport_matrix[i];                        rtmp4=rtmp*omega-fctp->main_diagonal_low_order_transport_matrix[i];
213                        if (ABS(rtmp3)>0) rtmp4+=sourceN[i]*rtmp2/rtmp3;                        if (ABS(rtmp3)>0) rtmp4+=sourceN[i]*rtmp2/rtmp3;
214                        fctp->iteration_matrix->mainBlock->val[fctp->main_iptr[i]]=rtmp4;                        fctp->iteration_matrix->mainBlock->val[fctp->main_iptr[i]]=rtmp4;
215                        uTilde_n[i]=rtmp3;                        uTilde_n[i]=rtmp3;
# Line 227  printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF Line 225  printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
225                        }                        }
226                        uTilde_n[i]=rtmp3;                        uTilde_n[i]=rtmp3;
227                   }                   }
228                     omega=1.;
229                   /* no update of iteration_matrix retquired! */                   /* no update of iteration_matrix retquired! */
230               } /* end (fctp->theta > 0) */               } /* end (fctp->theta > 0) */
231               /*               /*
# Line 241  printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF Line 240  printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
240                 */                 */
241                 m=0;                 m=0;
242                 converged=FALSE;                 converged=FALSE;
243                 while ( (!converged) && (m<50) && Paso_noError()) {                 while ( (!converged) && (m<500) && Paso_noError()) {
244                      printf("iteration step %d\n",m+1);                      printf("iteration step %d\n",m+1);
245                      /*                      /*
246                       *  set the ant diffusion fluxes:                       *  set the ant diffusion fluxes:
247                       *                       *
                      *   initially we set f_{ij} = - dt d_{ij} (u[j]-u[i])  
                      *   and then f_{ij} += omega (m_{ij} - dt (1-theta) d_{ij}) (du[j]-du[i])  
248                       */                       */
249                      if (m==0) {                       Paso_FCTransportProblem_setAntiDiffusionFlux(dt2,fctp,flux_matrix,u,fctp->u);
                        Paso_SystemMatrix_setValues(flux_matrix,0.);  
                        Paso_FCTransportProblem_updateAntiDiffusionFlux(fctp,flux_matrix,0.,-dt,u);  
                     } else {  
                        Paso_FCTransportProblem_updateAntiDiffusionFlux(fctp,flux_matrix,  
                                                                        omega,-omega*dt*(1.-fctp->theta),du_m);  
                     }  
250                      /*                      /*
251                       *  apply pre flux-correction: f_{ij}:=0 if f_{ij}*(\tilde{u}[i]- \tilde{u}[j])<=0                       *  apply pre flux-correction: f_{ij}:=0 if f_{ij}*(\tilde{u}[i]- \tilde{u}[j])<=0
252                       *                       *
253                       *  this is not entirely correct!!!!!                       *  this is not entirely correct!!!!!
254                       *                       *
255                       */                       */
256                      Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(flux_matrix,uTilde_n);                      Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(flux_matrix,uTilde_n);
257                      /*                      /*
258                       *  set flux limms RN_m,RP_m                       *  set flux limms RN_m,RP_m
259                       *                       *
260                       */                       */
261                      Paso_FCTransportProblem_setRs(flux_matrix,fctp->lumped_mass_matrix,QN_n,QP_n,RN_m,RP_m);                      Paso_FCTransportProblem_setRs(flux_matrix,fctp->lumped_mass_matrix,QN_n,QP_n,RN_m,RP_m);
262                      /*                      /*
263                       * z_m[i]=m_i*u[i] - dt*theta*sum_{j<>i} l_{ij} (u[j]-u[i]) + dt q^-[i]                       * z_m[i]=b_n[i] - (m_i*u[i] - dt2*theta*sum_{j<>i} l_{ij} (u[j]-u[i]) + dt2 q^-[i])
264                       *                       *
265                       * note that iteration_matrix stores the negative values of the                       * note that iteration_matrix stores the negative values of the
266                       * low order transport matrix l therefore a=dt*fctp->theta is used.                       * low order transport matrix l therefore a=dt2*fctp->theta is used.
267                       */                       */
268    
269                      Paso_SolverFCT_setMuPaLuPbQ(z_m,fctp->lumped_mass_matrix,u,                      Paso_SolverFCT_setMuPaLuPbQ(z_m,fctp->lumped_mass_matrix,u,
270                                                  dt*fctp->theta,fctp->iteration_matrix,dt,sourceN);                                                  dt2*fctp->theta,fctp->iteration_matrix,dt2,sourceN);
271                        #pragma omp parallel for private(i)
272                        for (i = 0; i < n_rows; ++i) z_m[i]=b_n[i]-z_m[i];
273    
274                       /* add corrected fluxes into z_m */                       /* add corrected fluxes into z_m */
275                       Paso_FCTransportProblem_addCorrectedFluxes(z_m,flux_matrix,RN_m,RP_m);                       Paso_FCTransportProblem_addCorrectedFluxes(z_m,flux_matrix,RN_m,RP_m);
276                       /*                       /*
277                        * now we solve the linear system to get the correction dt:                        * now we solve the linear system to get the correction dt2:
278                        *                        *
279                        */                        */
280                        if (fctp->theta > 0) {                        if (fctp->theta > 0) {
281                              /*  set the right hand side of the linear system: */                              /*  set the right hand side of the linear system: */
282                              #pragma omp parallel for private(i)                              options->tolerance=1.e-2;
                             for (i = 0; i < n_rows; ++i) {  
                                z_m[i]=b_n[i]-z_m[i];  
                             }  
                             options->tolerance=1.e-3;  
283                              Paso_solve(fctp->iteration_matrix,du_m,z_m,options);                              Paso_solve(fctp->iteration_matrix,du_m,z_m,options);
284                              /* TODO: check errors ! */                              /* TODO: check errors ! */
                             omega=dt*fctp->theta;  
285                        } else {                        } else {
286                              #pragma omp parallel for private(i,rtmp,rtmp1)                              #pragma omp parallel for private(i,rtmp,rtmp1)
287                              for (i = 0; i < n_rows; ++i) {                              for (i = 0; i < n_rows; ++i) {
288                                  rtmp=fctp->lumped_mass_matrix[i];                                  rtmp=fctp->lumped_mass_matrix[i];
289                                  if (ABS(rtmp)>0.) {                                  if (ABS(rtmp)>0.) {
290                                     rtmp1=(b_n[i]-z_m[i])/rtmp;                                     rtmp1=z_m[i]/rtmp;
291                                  } else {                                  } else {
292                                     rtmp1=0;                                     rtmp1=0;
293                                  }                                  }
294                                  du_m[i]=rtmp1;                                  du_m[i]=rtmp1;
295                              }                              }
                             omega=1.;  
296                        }                        }
297                        /*                        /*
298                         * update u and calculate norm of du_m and the new u:                         * update u and calculate norm of du_m and the new u:
# Line 340  printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF Line 328  printf("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
328                         m++;                         m++;
329                         printf("iteration step %d: norm u and du_m : %e %e\n",m,norm_u,norm_du);                         printf("iteration step %d: norm u and du_m : %e %e\n",m,norm_u,norm_du);
330                         /* TODO: check if du_m has been redu_mced */                         /* TODO: check if du_m has been redu_mced */
331                      } /* end of inner mation */                      } /* end of inner iteration */
332                        #pragma omp parallel for schedule(static) private(i)
333                        for (i=0;i<n_rows;++i) fctp->u[i]=u[i];
334                        n++;
335                 } /* end of time integration */                 } /* end of time integration */
336                 #pragma omp parallel for schedule(static) private(i)                 #pragma omp parallel for schedule(static) private(i)
337                 for (i=0;i<n_rows;++i) fctp->u[i]=u[i]+fctp->u_min;                 for (i=0;i<n_rows;++i) u[i]=fctp->u[i]+fctp->u_min;
338                 /* TODO: update u_min */                 /* TODO: update u_min ? */
339    
340          }          }
341          /*          /*

Legend:
Removed from v.1407  
changed lines
  Added in v.1410

  ViewVC Help
Powered by ViewVC 1.1.26