/[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 1415 by ksteube, Thu Feb 21 04:57:17 2008 UTC revision 1661 by gross, Mon Jul 21 22:08:27 2008 UTC
# Line 16  Line 16 
16  /* Paso: Flux correction transport solver  /* Paso: Flux correction transport solver
17   *   *
18   * solves Mu_t=Ku+q   * solves Mu_t=Ku+q
19   *        u(0) >= u_min   *        u(0) >=0
20   *   *
21   * Warning: the program assums sum_{j} k_{ij}=0!!!   * Warning: the program assums sum_{j} k_{ij}=0!!!
22   *   *
# Line 46  double Paso_FCTransportProblem_getSafeTi Line 46  double Paso_FCTransportProblem_getSafeTi
46     n_rows=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     n_rows=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
47     if (! fctp->valid_matrices) {     if (! fctp->valid_matrices) {
48          fctp->dt_max=LARGE_POSITIVE_FLOAT;          fctp->dt_max=LARGE_POSITIVE_FLOAT;
   
   
49          /* extract the row sum of the advective part */          /* extract the row sum of the advective part */
50          Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);          Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);
51    
# Line 65  double Paso_FCTransportProblem_getSafeTi Line 63  double Paso_FCTransportProblem_getSafeTi
63                 for (i=0;i<n_rows;++i) {                 for (i=0;i<n_rows;++i) {
64                      rtmp1=fctp->main_diagonal_low_order_transport_matrix[i];                      rtmp1=fctp->main_diagonal_low_order_transport_matrix[i];
65                      rtmp2=fctp->lumped_mass_matrix[i];                      rtmp2=fctp->lumped_mass_matrix[i];
66                      if ( (rtmp1<0 && rtmp2>=0.) || (rtmp1>0 && rtmp2<=0.) ) {                      if ( (rtmp1<0 && rtmp2>0.) || (rtmp1>0 && rtmp2<0) ) {
67                          dt_max_loc=MIN(dt_max_loc,-rtmp2/rtmp1);                          dt_max_loc=MIN(dt_max_loc,-rtmp2/rtmp1);
68                      }                      }
69                  }                  }
# Line 95  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;
97     index_t i, j;     long 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_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 113  void Paso_SolverFCT_solve(Paso_FCTranspo Line 114  void Paso_SolverFCT_solve(Paso_FCTranspo
114      *  allocate memory      *  allocate memory
115      *      *
116      */      */
117     Paso_SystemMatrix_allocBuffer(fctp->iteration_matrix);     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 134  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,
148                                         fctp->transport_matrix->col_block_size);                                         fctp->transport_matrix->col_block_size);
149     if (Paso_noError()) {     if (Paso_noError()) {
        Paso_SystemMatrix_allocBuffer(flux_matrix);  
        Paso_SystemMatrix_allocBuffer(fctp->iteration_matrix);  
150         /*         /*
151          *    Preparation:          *    Preparation:
152          *          *
# Line 148  void Paso_SolverFCT_solve(Paso_FCTranspo Line 154  void Paso_SolverFCT_solve(Paso_FCTranspo
154            
155         /* decide on substepping */         /* decide on substepping */
156         if (fctp->dt_max < LARGE_POSITIVE_FLOAT) {         if (fctp->dt_max < LARGE_POSITIVE_FLOAT) {
157            n_substeps=ceil(dt/dt_max);            n_substeps=(long)ceil(dt/dt_max);
158         } else {         } else {
159            n_substeps=1;            n_substeps=1;
160         }         }
161         dt2=dt/n_substeps;         dt2=dt/n_substeps;
162         printf("%d time steps of size is %e (theta = %e, dt_max=%e).\n",n_substeps, dt2,fctp->theta, dt_max);         printf("%ld time steps of size is %e (theta = %e, dt_max=%e).\n",n_substeps, dt2,fctp->theta, dt_max);
163         /*         /*
164      * seperate source into positive and negative part:      * seperate source into positive and negative part:
165      */      */
# Line 168  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 177  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              printf("substep step %d at t=%e\n",n+1,t);  
188              #pragma omp parallel for private(i)              printf("substep step %ld at t=%e\n",n+1,t);
              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 189  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 207  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 221  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 %ld\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 266  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 306  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 316  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                  double local_norm[2], norm[2];
343                            local_norm[0]=norm_u;                            local_norm[0]=norm_u;
344                            local_norm[1]=norm_du;                            local_norm[1]=norm_du;
345                        MPI_Allreduce(local_norm,norm, 2, MPI_DOUBLE, MPI_MAX, fctp->mpi_info->comm);                        MPI_Allreduce(local_norm,norm, 2, MPI_DOUBLE, MPI_MAX, fctp->mpi_info->comm);
# Line 326  void Paso_SolverFCT_solve(Paso_FCTranspo Line 349  void Paso_SolverFCT_solve(Paso_FCTranspo
349                         norm_du*=omega;                         norm_du*=omega;
350                         converged=(norm_du <= tolerance * norm_u);                         converged=(norm_du <= tolerance * norm_u);
351                         m++;                         m++;
352                         printf("iteration step %d: norm u and du_m : %e %e\n",m,norm_u,norm_du);                         printf("iteration step %ld: norm u and du_m : %e %e\n",m,norm_u,norm_du);
353                         /* TODO: check if du_m has been redu_mced */                         /* TODO: check if du_m has been redu_mced */
354                           Paso_Coupler_finishCollect(u_m_coupler);
355                      } /* end of inner iteration */                      } /* end of inner iteration */
356                      #pragma omp parallel for schedule(static) private(i)                      SWAP(fctp->u,u_m,double*);
357                      for (i=0;i<n_rows;++i) fctp->u[i]=u[i];                      SWAP(fctp->u_coupler,u_m_coupler,Paso_Coupler*);
358                      n++;                      n++;
359                 } /* end of time integration */                 } /* end of time integration */
360                   options->tolerance=tolerance;
361                 #pragma omp parallel for schedule(static) private(i)                 #pragma omp parallel for schedule(static) private(i)
362                 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];
                /* TODO: update u_min ? */  
   
363          }          }
364          /*          /*
365           *  clean-up:           *  clean-up:
# Line 353  void Paso_SolverFCT_solve(Paso_FCTranspo Line 376  void Paso_SolverFCT_solve(Paso_FCTranspo
376          TMPMEMFREE(RP_m);          TMPMEMFREE(RP_m);
377          TMPMEMFREE(z_m);          TMPMEMFREE(z_m);
378          TMPMEMFREE(du_m);          TMPMEMFREE(du_m);
379            TMPMEMFREE(u_m);
380            Paso_Coupler_free(QN_n_coupler);
381            Paso_Coupler_free(QP_n_coupler);
382            Paso_Coupler_free(RN_m_coupler);
383            Paso_Coupler_free(RP_m_coupler);
384            Paso_Coupler_free(uTilde_n_coupler);
385            Paso_Coupler_free(u_m_coupler);
386  }  }

Legend:
Removed from v.1415  
changed lines
  Added in v.1661

  ViewVC Help
Powered by ViewVC 1.1.26