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

revision 1370 by gross, Wed Jan 2 09:21:43 2008 UTC revision 1400 by gross, Thu Jan 24 06:04:31 2008 UTC
# Line 62  Line 62
62
63          For the case theta==0 (explicit case) no iteration is required. One sets          For the case theta==0 (explicit case) no iteration is required. One sets
64
65                M*u= M*u_last + dt*F(u_last)                M*u= M*u_last + dt*b
66
67          with F(u) =  L u + f_a(u).          with b=F(u) =  L u + f_a(u).
68
69
70          For the case theta>0 we set          For the case theta>0 we solve (*) is the form
71
72                [L-M/theta*dt] du2 =M/(theta*dt)*u_last + ((1-theta)/theta)*F(u_last) - M/(theta*dt)*u + F(u)
73
74            for du2=-du which is solved as
75
76                A du2  = r= b + f(u)
77
78            with
79
80                du=-du2
81              A=L-M/(theta*dt) (stored into transport_matrix)              A=L-M/(theta*dt) (stored into transport_matrix)
82              F(u)=L u + f_a(u) = M/(theta*dt) u + A u + f_a(u)              f(u) =-M/(theta*dt)*u + F(u) = - M/(theta*dt)*u + Lu + f_a(u) = f_a(u) + Au
83              b(u)= - M/(theta*dt) * u + F(u) =  A u + f_a(u)              F(u)=f(u)+M/(theta*dt)*u
84              b_last=M/(theta*dt)*u_last + (1-theta)*dt * F(u_last)              b= M/(theta*dt)*u_last+(1-theta)/theta*F(u_last)=
85                    =M/(theta*dt)*u_last + (1-theta)/theta * [ M/(theta*dt) u_last + A u_last + f_a(u_last) ]               = M/(theta*dt)*u_last+(1-theta)/theta*(f(u)+M/(theta*dt)*u_last)
86                    =(1/theta**2*dt)*M u_last + (1-theta)/theta*b(u_last)               = M/(theta*dt)*u_last+(1-theta)/theta*(f(u_last)+M/(theta*dt)*u_last)
87          so (*) takes the form               = M*1./(theta*dt)*(1+(1-theta)/theta)*u_last+(1-theta)/theta*f(u_last)

b_last=(1/theta**2*dt)*M u_last + (1-theta)/theta*b(u_last)
while (|du| > tol * |u|) :
A*du=b_last + b(u)
u <- u-du
88
89  */  */
90
91  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) {
92
93     index_t i;     index_t i, j;
94     int n_substeps,n;     int n_substeps,n, iter;
95       double fac, fac2, *b=NULL, *f=NULL, *du=NULL;
96     double dt2=fctp->dt_max, dt2_loc, rtmp,rtmp2,t;     double dt2=fctp->dt_max, dt2_loc, rtmp,rtmp2,t;
97     printf("FFFF %e",fctp->dt_max);     double local_norm[2],norm[2],local_norm_u,local_norm_du,norm_u,norm_du, tolerance;
98       dim_t n_rows=Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);
99       bool_t converged;
100     if (dt<=0.) {     if (dt<=0.) {
101         Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");         Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");
102     }     }
# Line 101  void Paso_SolverFCT_solve(Paso_FCTranspo Line 108  void Paso_SolverFCT_solve(Paso_FCTranspo
109          /* create a copy of the main diagonal entires of the transport-matrix */          /* create a copy of the main diagonal entires of the transport-matrix */
110          #pragma omp parallel for schedule(static) private(i)          #pragma omp parallel for schedule(static) private(i)
111          for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {          for (i=0;i<n_rows;++i) {
112                 fctp->transport_matrix_diagonal[i]=                 fctp->transport_matrix_diagonal[i]=
113                          fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];                          fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];
114          }          }
115
116  Paso_SystemMatrix_saveMM(fctp->flux_matrix,"flux.mm");   /* Paso_SystemMatrix_saveMM(fctp->flux_matrix,"flux.mm");
117  Paso_SystemMatrix_saveMM(fctp->transport_matrix,"trans.mm");   Paso_SystemMatrix_saveMM(fctp->transport_matrix,"trans.mm");  */
118          /*=================================================================== *          /*=================================================================== *
119           *           *
120           *  calculate time step size:                                                     *  calculate time step size:
121           */           */
122          dt2=fctp->dt_max;          dt2=LARGE_POSITIVE_FLOAT;
123          if (fctp->theta < 1.) {          if (fctp->theta < 1.) {
dt2=LARGE_POSITIVE_FLOAT;
124              #pragma omp parallel private(dt2_loc)              #pragma omp parallel private(dt2_loc)
125              {              {
126                 dt2_loc=LARGE_POSITIVE_FLOAT;                 dt2_loc=LARGE_POSITIVE_FLOAT;
127                 #pragma omp for schedule(static) private(i,rtmp,rtmp2)                 #pragma omp for schedule(static) private(i,rtmp,rtmp2)
128                 for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {                 for (i=0;i<n_rows;++i) {
129                      rtmp=fctp->transport_matrix_diagonal[i];                      rtmp=fctp->transport_matrix_diagonal[i];
130                      rtmp2=fctp->lumped_mass_matrix[i];                      rtmp2=fctp->lumped_mass_matrix[i];
131                      if ( (rtmp<0 && rtmp2>=0.) || (rtmp>0 && rtmp2<=0.) ) {                      if ( (rtmp<0 && rtmp2>=0.) || (rtmp>0 && rtmp2<=0.) ) {
132                          dt2_loc=MIN(dt2_loc,-rtmp2/rtmp);                          dt2_loc=MIN(dt2_loc,-rtmp2/rtmp);
133                      }                      }
134                  }                  }
135                  #pragma omp critcal                  #pragma omp critical
136                  {                  {
137                      dt2=MIN(dt2,dt2_loc);                      dt2=MIN(dt2,dt2_loc);
138                  }                  }
# Line 135  Paso_SystemMatrix_saveMM(fctp->transport Line 141  Paso_SystemMatrix_saveMM(fctp->transport
141                 dt2_loc = dt2;                 dt2_loc = dt2;
142             MPI_Allreduce(&dt2_loc, &dt2, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);             MPI_Allreduce(&dt2_loc, &dt2, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
143              #endif              #endif
144              dt2*=1./(1.-fctp->theta);              if (dt2<LARGE_POSITIVE_FLOAT) dt2*=1./(1.-fctp->theta);
if (fctp->dt_max>0.) dt2=MIN(dt2,fctp->dt_max);
}
if (dt2 > 0.) {
fctp->dt=dt2;
} else {
fctp->dt=fctp->dt_max;
145          }          }
146          if (fctp->dt < 0.) {          if (dt2 <= 0.) {
147              Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt_max must be positive.");              Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");
148          } else {          } else {
149             printf("minimum time step size is %e (theta = %e).\n",fctp->dt,fctp->theta);               if (dt2<LARGE_POSITIVE_FLOAT) printf("minimum time step size is %e (theta = %e).\n",dt2,fctp->theta);
}
/* ===========================================================
*
*
*/
if (Paso_noError()) {

150          }          }
151            fctp->dt=dt2;
152            fctp->dt_max=dt2; /* FIXME: remove*/
153          fctp->valid_matrices=Paso_noError();          fctp->valid_matrices=Paso_noError();
154     }     }
155
/* b_last=M*u_last + (1-theta) * F(u_last) */

/* decide on substepping */
n_substeps=ceil(dt/fctp->dt);
dt2=dt/n_substeps;
printf("%d time steps of size is %e (theta = %e, dt_max=%e).\n",n_substeps, dt2,fctp->theta,fctp->dt);

156     /* */     /* */
157     Paso_SystemMatrix_allocBuffer(fctp->transport_matrix);     Paso_SystemMatrix_allocBuffer(fctp->transport_matrix);
158       /*
159        * allocate memory:
160        *
161        */
162       b=MEMALLOC(n_rows,double);
163       Paso_checkPtr(b);
164     if (fctp->theta>0) {     if (fctp->theta>0) {
165          #pragma omp parallel for schedule(static) private(i)         b=MEMALLOC(n_rows,double);
166          for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {         du=MEMALLOC(n_rows,double);
167                 fctp->transport_matrix_diagonal[i]=         f=MEMALLOC(n_rows,double);
168                          fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];         Paso_checkPtr(du);
169          }         Paso_checkPtr(f);
170     } else {     }
171        /* u= u_last + M^-1*dt*F(u_last) */     if (Paso_noError()) {
172        t=dt2;         /*
173        n=0;          *    Preparation:
174        while(n<n_substeps) {          *
175          printf("substep step %d at t=%e\n",n+1,t);          */
176          Paso_FCTransportProblem_setFlux(fctp,fctp->u,u); /* u stores F(u_last)*/
177          #pragma omp parallel for schedule(static) private(i)         /* decide on substepping */
178          for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {         if (fctp->dt < LARGE_POSITIVE_FLOAT) {
179              rtmp=fctp->u[i];            n_substeps=ceil(dt/fctp->dt);
180              fctp->u[i]=fctp->u[i]+dt2*u[i]/fctp->lumped_mass_matrix[i];         } else {
181              printf("%d : %e %e %e\n",i,u[i],fctp->u[i],rtmp);            n_substeps=1.;
182           }
183           dt2=dt/n_substeps;
184           printf("%d time steps of size is %e (theta = %e, dt_max=%e).\n",n_substeps, dt2,fctp->theta,fctp->dt);
185
186           /*
187            *  implicit case:
188            *
189            *   A=L-M/(theta*dt) (stored into transport_matrix)
190            *
191            * b= M/(theta*dt)*u_last+(1-theta)/theta)*F(u_last)
192            *
193            */
194           if (fctp->theta>0) {
195              Paso_solve_free(fctp->transport_matrix);
196              fac=1./(fctp->theta*dt2);
197              fac2=(1.-fctp->theta)/fctp->theta;
198              #pragma omp parallel for schedule(static) private(i)
199              for (i=0;i<n_rows;++i) {
200                    fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]]=
201                                                fctp->transport_matrix_diagonal[i]-fctp->lumped_mass_matrix[i]*fac;
202              }
203           }
204           #pragma omp parallel for schedule(static) private(i)
205           for (i=0;i<n_rows;++i) {
206             u[i]=fctp->u[i];
207    /* printf("A %d : %e %e\n",i,u[i],fctp->u[i]); */
208           }
209           /*
210            * now the show can begin:
211            *
212            */
213            t=dt2;
214            n=0;
215            tolerance=options->tolerance;
216            while(n<n_substeps && Paso_noError()) {
217                printf("substep step %d at t=%e\n",n+1,t);
218                if (fctp->theta>0.) {
219                    /*
220                     * implicit scheme:
221                     *
222                     */
223                    if (fctp->theta<1.) {
224                       Paso_FCTransportProblem_setFlux(fctp,u,b);  /* b=f(u_last) */
225                       #pragma omp parallel for schedule(static) private(i)
226                       for (i=0;i<n_rows;++i) {
227                            b[i]=fctp->lumped_mass_matrix[i]*fac*(1.+fac2)*u[i]+fac2*b[i];
228                        }
229                    } else {
230                       #pragma omp parallel for schedule(static) private(i)
231                       for (i=0;i<n_rows;++i) {
232                           b[i]=fctp->lumped_mass_matrix[i]*fac*u[i];
233                       }
234                    }
235                    /*
236                     * Enter iteration on a time step :
237                     *
238                     */
239                     iter=0;
240                     converged=FALSE;
241                     while ( (!converged) && (iter<50) && Paso_noError()) {
242                        printf("iteration step %d\n",iter+1);
243                        Paso_FCTransportProblem_setFlux(fctp,u,f);
244                        #pragma omp parallel for schedule(static) private(i)
245                        for (i=0;i<n_rows;++i) {
246                           f[i]+=b[i];
247                        }
248                        options->tolerance=1.e-3;
249                        Paso_solve(fctp->transport_matrix,du,f,options);
250                        /* update u and calculate norms */
251                        norm_u=0.;
252                        norm_du=0.;
253                        #pragma omp parallel private(local_norm_u,local_norm_du)
254                        {
255                            local_norm_u=0.;
256                            local_norm_du=0.;
257                            #pragma omp for schedule(static) private(i)
258                            for (i=0;i<n_rows;++i) {
259                                 u[i]-=du[i];
260                                 local_norm_u=MAX(local_norm_u,ABS(u[i]));
261                                 local_norm_du=MAX(local_norm_du,ABS(du[i]));
262                            }
263                            #pragma omp crtical
264                            {
265                                norm_u=MAX(norm_u,local_norm_u);
266                                norm_du=MAX(norm_du,local_norm_du);
267                            }
268                        }
269                        #ifdef PASO_MPI
270                           local_norm[0]=norm_u;
271                           local_norm[1]=norm_du;
272                       MPI_Allreduce(local_norm,norm, 2, MPI_DOUBLE, MPI_MAX, fctp->mpi_info->comm);
273                           norm_u=norm[0];
274                           norm_du=norm[1];
275                        #endif
276                        converged=(norm_du <= tolerance * norm_u);
277                        iter++;
278                        printf("iteration step %d: norm u and du : %e %e\n",iter,norm_u,norm_du);
279                     }
280                 } else {
281                    /*
282                      * explicit scheme:
283                      *
284                     */
285                    #pragma omp parallel for schedule(static) private(i)
286                    for (i=0;i<n_rows;++i) {
287                        u[i]+=dt2*b[i]/fctp->lumped_mass_matrix[i];
288                    }
289                 }
290                 /* and the next time step */
291                 t+=dt2;
292                 n++;
293          }          }
294          t+=dt2;          /*
295          n++;           * save last u
296        }           *
297    }           */
298    if (Paso_noError()) {          if (Paso_noError()) {
299       #pragma omp parallel for schedule(static) private(i)               #pragma omp parallel for schedule(static) private(i)
300       for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {               for (i=0;i<n_rows;++i) {
301         u[i]=fctp->u[i];                   fctp->u[i]=u[i];
302       }               }
303    }          }
304  }            }
305        /*
306         *
307         * clean-up
308         *
309         */
310        if (fctp->theta>0) {
311                MEMFREE(b);
312                MEMFREE(du);
313                MEMFREE(f);
314        }
315    }

Legend:
 Removed from v.1370 changed lines Added in v.1400