/[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 1561 by gross, Thu May 8 08:52:41 2008 UTC revision 1562 by gross, Wed May 21 13:04:40 2008 UTC
# Line 93  double Paso_FCTransportProblem_getSafeTi Line 93  double Paso_FCTransportProblem_getSafeTi
93    
94    
95  void Paso_SolverFCT_solve(Paso_FCTransportProblem* fctp, double* u, double dt, double* source, Paso_Options* options) {  void Paso_SolverFCT_solve(Paso_FCTransportProblem* fctp, double* u, double dt, double* source, Paso_Options* options) {
   
96     index_t i, j;     index_t i, j;
97     int n_substeps,n, m;     int n_substeps,n, m;
98     double dt_max, omega, dt2,t;     double dt_max, omega, dt2,t;
99     double local_norm[2],norm[2],local_norm_u,local_norm_du,norm_u,norm_du, tolerance;     double local_norm[2],norm[2],local_norm_u,local_norm_du,norm_u,norm_du, tolerance;
100     register double rtmp1,rtmp2,rtmp3,rtmp4, rtmp;     register double rtmp1,rtmp2,rtmp3,rtmp4, rtmp;
101     double *b_n=NULL, *sourceP=NULL, *sourceN=NULL, *uTilde_n=NULL, *QN_n=NULL, *QP_n=NULL, *RN_m=NULL, *RP_m=NULL, *z_m=NULL, *du_m=NULL;     double *b_n=NULL, *sourceP=NULL, *sourceN=NULL, *uTilde_n=NULL, *QN_n=NULL, *QP_n=NULL, *RN_m=NULL, *RP_m=NULL, *z_m=NULL, *du_m=NULL,*u_m=NULL;
102       Paso_Coupler* QN_n_coupler=NULL, *QP_n_coupler=NULL, *uTilde_n_coupler=NULL, *RN_m_coupler=NULL, *RP_m_coupler=NULL, *u_m_coupler=NULL;
103     Paso_SystemMatrix *flux_matrix=NULL;     Paso_SystemMatrix *flux_matrix=NULL;
104     dim_t n_rows=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     dim_t n_rows=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
105     bool_t converged;     bool_t converged;
106       dim_t max_m=500;
107       double inner_tol=1.e-2;
108       dim_t blockSize=Paso_FCTransportProblem_getBlockSize(fctp);
109     if (dt<=0.) {     if (dt<=0.) {
110         Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");         Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");
111     }     }
# Line 111  void Paso_SolverFCT_solve(Paso_FCTranspo Line 114  void Paso_SolverFCT_solve(Paso_FCTranspo
114      *  allocate memory      *  allocate memory
115      *      *
116      */      */
117       u_m=TMPMEMALLOC(n_rows,double);
118       Paso_checkPtr(u_m);
119     b_n=TMPMEMALLOC(n_rows,double);     b_n=TMPMEMALLOC(n_rows,double);
120     Paso_checkPtr(b_n);     Paso_checkPtr(b_n);
121     sourceP=TMPMEMALLOC(n_rows,double);     sourceP=TMPMEMALLOC(n_rows,double);
# Line 131  void Paso_SolverFCT_solve(Paso_FCTranspo Line 136  void Paso_SolverFCT_solve(Paso_FCTranspo
136     Paso_checkPtr(z_m);     Paso_checkPtr(z_m);
137     du_m=TMPMEMALLOC(n_rows,double);     du_m=TMPMEMALLOC(n_rows,double);
138     Paso_checkPtr(du_m);     Paso_checkPtr(du_m);
139       QN_n_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(fctp),blockSize);
140       QP_n_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(fctp),blockSize);
141       RN_m_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(fctp),blockSize);
142       RP_m_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(fctp),blockSize);
143       uTilde_n_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(fctp),blockSize);
144       u_m_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(fctp),blockSize);
145     flux_matrix=Paso_SystemMatrix_alloc(fctp->transport_matrix->type,     flux_matrix=Paso_SystemMatrix_alloc(fctp->transport_matrix->type,
146                                         fctp->transport_matrix->pattern,                                         fctp->transport_matrix->pattern,
147                                         fctp->transport_matrix->row_block_size,                                         fctp->transport_matrix->row_block_size,
# Line 163  void Paso_SolverFCT_solve(Paso_FCTranspo Line 174  void Paso_SolverFCT_solve(Paso_FCTranspo
174               sourceP[i]= rtmp;               sourceP[i]= rtmp;
175            }            }
176          }          }
177            u_m_coupler->data=u_m;
178  /*        for (i = 0; i < n_rows; ++i) printf("%d : %e \n",i,source[i],sourceP[i],sourceN[i]) */  /*        for (i = 0; i < n_rows; ++i) printf("%d : %e \n",i,source[i],sourceP[i],sourceN[i]) */
179          /*          /*
180           * let the show begin!!!!           * let the show begin!!!!
# Line 172  void Paso_SolverFCT_solve(Paso_FCTranspo Line 184  void Paso_SolverFCT_solve(Paso_FCTranspo
184          n=0;          n=0;
185          tolerance=options->tolerance;          tolerance=options->tolerance;
186          while(n<n_substeps && Paso_noError()) {          while(n<n_substeps && Paso_noError()) {
187    
188              printf("substep step %d at t=%e\n",n+1,t);              printf("substep step %d at t=%e\n",n+1,t);
             #pragma omp parallel for private(i)  
              for (i = 0; i < n_rows; ++i) {  
                       u[i]=fctp->u[i];  
              }  
189              /*              /*
190               * b^n[i]=m u^n[i] + dt2*(1-theta) sum_{j <> i} l_{ij}*(u^n[j]-u^n[i]) + dt2*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]
191               *               *
# Line 184  void Paso_SolverFCT_solve(Paso_FCTranspo Line 193  void Paso_SolverFCT_solve(Paso_FCTranspo
193               * low order transport matrix l therefore a=-dt2*(1-fctp->theta) is used.               * low order transport matrix l therefore a=-dt2*(1-fctp->theta) is used.
194               *               *
195               */               */
196               Paso_SolverFCT_setMuPaLuPbQ(b_n,fctp->lumped_mass_matrix,u,               Paso_SolverFCT_setMuPaLuPbQ(b_n,fctp->lumped_mass_matrix,fctp->u_coupler,
197                                            -dt2*(1-fctp->theta),fctp->iteration_matrix,dt2,sourceP);                                            -dt2*(1-fctp->theta),fctp->iteration_matrix,dt2,sourceP);
198               /*               /*
199            *   uTilde_n[i]=b[i]/m[i]            *   uTilde_n[i]=b[i]/m[i]
200            *            *
201            *   a[i,i]=m[i]/(dt2 theta) + \frac{1}{\theta} \frac{q^-[i]}-l[i,i]            *   fctp->iteration_matrix[i,i]=m[i]/(dt2 theta) + \frac{1}{\theta} \frac{q^-[i]}-l[i,i]
202            *            *
203            */            */
204               if (fctp->theta > 0) {               if (fctp->theta > 0) {
# Line 202  void Paso_SolverFCT_solve(Paso_FCTranspo Line 211  void Paso_SolverFCT_solve(Paso_FCTranspo
211                        if (ABS(rtmp)>0.) {                        if (ABS(rtmp)>0.) {
212                           rtmp3=b_n[i]/rtmp;                           rtmp3=b_n[i]/rtmp;
213                        } else {                        } else {
214                           rtmp3=u[i];                           rtmp3=fctp->u[i];
215                        }                        }
216                        rtmp4=rtmp*omega-fctp->main_diagonal_low_order_transport_matrix[i];                        rtmp4=rtmp*omega-fctp->main_diagonal_low_order_transport_matrix[i];
217                        if (ABS(rtmp3)>0) rtmp4+=sourceN[i]*rtmp2/rtmp3;                        if (ABS(rtmp3)>0) rtmp4+=sourceN[i]*rtmp2/rtmp3;
# Line 216  void Paso_SolverFCT_solve(Paso_FCTranspo Line 225  void Paso_SolverFCT_solve(Paso_FCTranspo
225                        if (ABS(rtmp)>0.) {                        if (ABS(rtmp)>0.) {
226                           rtmp3=b_n[i]/rtmp;                           rtmp3=b_n[i]/rtmp;
227                        } else {                        } else {
228                           rtmp3=u[i];                           rtmp3=fctp->u[i];
229                        }                        }
230                        uTilde_n[i]=rtmp3;                        uTilde_n[i]=rtmp3;
231                   }                   }
232                   omega=1.;                   omega=1.;
233                   /* no update of iteration_matrix retquired! */                   /* no update of iteration_matrix required! */
234               } /* end (fctp->theta > 0) */               } /* end (fctp->theta > 0) */
235    
236                 /* distribute uTilde_n: */
237                 Paso_Coupler_startCollect(uTilde_n_coupler,uTilde_n);
238                 Paso_Coupler_finishCollect(uTilde_n_coupler);
239               /*               /*
240                * calculate QP_n[i] max_{j} (\tilde{u}[j]- \tilde{u}[i] )                * calculate QP_n[i] max_{j} (\tilde{u}[j]- \tilde{u}[i] )
241                *           QN_n[i] min_{j} (\tilde{u}[j]- \tilde{u}[i] )                *           QN_n[i] min_{j} (\tilde{u}[j]- \tilde{u}[i] )
242                *                *
243                */                */
244                Paso_SolverFCT_setQs(uTilde_n,QN_n,QP_n,fctp->iteration_matrix);                Paso_SolverFCT_setQs(uTilde_n_coupler,QN_n,QP_n,fctp->iteration_matrix);
245                  Paso_Coupler_startCollect(QN_n_coupler,QN_n);
246                  Paso_Coupler_startCollect(QP_n_coupler,QP_n);
247                  Paso_Coupler_finishCollect(QN_n_coupler);
248                  Paso_Coupler_finishCollect(QP_n_coupler);
249                /*                /*
250                 * now we enter the mation on a time-step:                 * now we enter the iteration on a time-step:
251                 *                 *
252                 */                 */
253                   Paso_Coupler_copyAll(fctp->u_coupler, u_m_coupler);
254    /* REPLACE BY JACOBEAN FREE NEWTON */
255                 m=0;                 m=0;
256                 converged=FALSE;                 converged=FALSE;
257                 while ( (!converged) && (m<500) && Paso_noError()) {                 while ( (!converged) && (m<max_m) && Paso_noError()) {
258                      printf("iteration step %d\n",m+1);                      printf("iteration step %d\n",m+1);
259                      /*                      /*
260                       *  set the ant diffusion fluxes:                       *  set the ant diffusion fluxes:
261                       *                       *
262                       */                      */
263                       Paso_FCTransportProblem_setAntiDiffusionFlux(dt2,fctp,flux_matrix,u,fctp->u);                      Paso_FCTransportProblem_setAntiDiffusionFlux(dt2,fctp,flux_matrix,u_m_coupler);
264                      /*                      /*
265                       *  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
266                       *                       *
                      *  this is not entirely correct!!!!!  
                      *  
267                       */                       */
268                      Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(flux_matrix,uTilde_n);                      Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(flux_matrix,uTilde_n_coupler);
269                      /*                      /*
270                       *  set flux limms RN_m,RP_m                       *  set flux limiters RN_m,RP_m
271                       *                       *
272                       */                       */
273                      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_coupler,QP_n_coupler,RN_m,RP_m);
274                        Paso_Coupler_startCollect(RN_m_coupler,RN_m);
275                        Paso_Coupler_startCollect(RP_m_coupler,RP_m);
276                      /*                      /*
277                       * z_m[i]=b_n[i] - (m_i*u[i] - dt2*theta*sum_{j<>i} l_{ij} (u[j]-u[i]) + dt2 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])
278                       *                       *
# Line 261  void Paso_SolverFCT_solve(Paso_FCTranspo Line 280  void Paso_SolverFCT_solve(Paso_FCTranspo
280                       * low order transport matrix l therefore a=dt2*fctp->theta is used.                       * low order transport matrix l therefore a=dt2*fctp->theta is used.
281                       */                       */
282    
283                      Paso_SolverFCT_setMuPaLuPbQ(z_m,fctp->lumped_mass_matrix,u,                      Paso_SolverFCT_setMuPaLuPbQ(z_m,fctp->lumped_mass_matrix, u_m_coupler,
284                                                  dt2*fctp->theta,fctp->iteration_matrix,dt2,sourceN);                                                  dt2*fctp->theta,fctp->iteration_matrix,dt2,sourceN);
285    /* for (i = 0; i < n_rows; ++i) {
286    * printf("%d: %e %e\n",i,z_m[i], b_n[i]);
287    * }
288    */
289                      #pragma omp parallel for private(i)                      #pragma omp parallel for private(i)
290                      for (i = 0; i < n_rows; ++i) z_m[i]=b_n[i]-z_m[i];                      for (i = 0; i < n_rows; ++i) z_m[i]=b_n[i]-z_m[i];
291    
292                       /* add corrected fluxes into z_m */                      Paso_Coupler_finishCollect(RN_m_coupler);
293                       Paso_FCTransportProblem_addCorrectedFluxes(z_m,flux_matrix,RN_m,RP_m);                      Paso_Coupler_finishCollect(RP_m_coupler);
294                       /*                      /* add corrected fluxes into z_m */
295                        Paso_FCTransportProblem_addCorrectedFluxes(z_m,flux_matrix,RN_m_coupler,RP_m_coupler);
296                        /*
297                        * now we solve the linear system to get the correction dt2:                        * now we solve the linear system to get the correction dt2:
298                        *                        *
299                        */                       */
300                        if (fctp->theta > 0) {                      if (fctp->theta > 0) {
301                              /*  set the right hand side of the linear system: */                              /*  set the right hand side of the linear system: */
302                              options->tolerance=1.e-2;                              options->tolerance=inner_tol;
303                                options->verbose=TRUE;
304                              Paso_solve(fctp->iteration_matrix,du_m,z_m,options);                              Paso_solve(fctp->iteration_matrix,du_m,z_m,options);
305                              /* TODO: check errors ! */                              /* TODO: check errors ! */
306                        } else {                        } else {
# Line 301  void Paso_SolverFCT_solve(Paso_FCTranspo Line 327  void Paso_SolverFCT_solve(Paso_FCTranspo
327                             local_norm_du=0.;                             local_norm_du=0.;
328                             #pragma omp for schedule(static) private(i)                             #pragma omp for schedule(static) private(i)
329                             for (i=0;i<n_rows;++i) {                             for (i=0;i<n_rows;++i) {
330                                  u[i]+=omega*du_m[i];                                  u_m[i]+=omega*du_m[i];
331                                  local_norm_u=MAX(local_norm_u,ABS(u[i]));                                  local_norm_u=MAX(local_norm_u,ABS(u_m[i]));
332                                  local_norm_du=MAX(local_norm_du,ABS(du_m[i]));                                  local_norm_du=MAX(local_norm_du,ABS(du_m[i]));
333                             }                             }
334                             #pragma omp critical                             #pragma omp critical
# Line 311  void Paso_SolverFCT_solve(Paso_FCTranspo Line 337  void Paso_SolverFCT_solve(Paso_FCTranspo
337                                 norm_du=MAX(norm_du,local_norm_du);                                 norm_du=MAX(norm_du,local_norm_du);
338                             }                             }
339                         }                         }
340                           Paso_Coupler_startCollect(u_m_coupler,u_m);
341                         #ifdef PASO_MPI                         #ifdef PASO_MPI
342                            local_norm[0]=norm_u;                            local_norm[0]=norm_u;
343                            local_norm[1]=norm_du;                            local_norm[1]=norm_du;
# Line 323  void Paso_SolverFCT_solve(Paso_FCTranspo Line 350  void Paso_SolverFCT_solve(Paso_FCTranspo
350                         m++;                         m++;
351                         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);
352                         /* TODO: check if du_m has been redu_mced */                         /* TODO: check if du_m has been redu_mced */
353                           Paso_Coupler_finishCollect(u_m_coupler);
354                      } /* end of inner iteration */                      } /* end of inner iteration */
355                      #pragma omp parallel for schedule(static) private(i)                      SWAP(fctp->u,u_m,double*);
356                      for (i=0;i<n_rows;++i) fctp->u[i]=u[i];                      SWAP(fctp->u_coupler,u_m_coupler,Paso_Coupler*);
357                      n++;                      n++;
358                 } /* end of time integration */                 } /* end of time integration */
359                   options->tolerance=tolerance;
360                 #pragma omp parallel for schedule(static) private(i)                 #pragma omp parallel for schedule(static) private(i)
361                 for (i=0;i<n_rows;++i) u[i]=fctp->u[i]+fctp->u_min;                 for (i=0;i<n_rows;++i) u[i]=fctp->u[i]+fctp->u_min;
                /* TODO: update u_min ? */  
   
362          }          }
363          /*          /*
364           *  clean-up:           *  clean-up:
# Line 348  void Paso_SolverFCT_solve(Paso_FCTranspo Line 375  void Paso_SolverFCT_solve(Paso_FCTranspo
375          TMPMEMFREE(RP_m);          TMPMEMFREE(RP_m);
376          TMPMEMFREE(z_m);          TMPMEMFREE(z_m);
377          TMPMEMFREE(du_m);          TMPMEMFREE(du_m);
378            TMPMEMFREE(u_m);
379            Paso_Coupler_free(QN_n_coupler);
380            Paso_Coupler_free(QP_n_coupler);
381            Paso_Coupler_free(RN_m_coupler);
382            Paso_Coupler_free(RP_m_coupler);
383            Paso_Coupler_free(uTilde_n_coupler);
384            Paso_Coupler_free(u_m_coupler);
385  }  }

Legend:
Removed from v.1561  
changed lines
  Added in v.1562

  ViewVC Help
Powered by ViewVC 1.1.26