/[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 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC revision 1401 by gross, Fri Jan 25 04:31:18 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, j;     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       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);     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 106  void Paso_SolverFCT_solve(Paso_FCTranspo Line 113  void Paso_SolverFCT_solve(Paso_FCTranspo
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;
# 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<n_rows;++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        /* n_substeps=1; */     if (Paso_noError()) {
172        /* u= u_last + M^-1*dt*F(u_last) */         /*
173        t=dt2;          *    Preparation:
174        n=0;          *
175        while(n<n_substeps) {          */
176          printf("substep step %d at t=%e\n",n+1,t);      
177          Paso_FCTransportProblem_setFlux(fctp,fctp->u,u); /* u stores F(u_last)*/         /* decide on substepping */
178          #pragma omp parallel for schedule(static) private(i)         if (fctp->dt < LARGE_POSITIVE_FLOAT) {
179          for (i=0;i<n_rows;++i) {            n_substeps=ceil(dt/fctp->dt);
180              rtmp=fctp->u[i];         } else {
181              fctp->u[i]=fctp->u[i]+dt2*u[i]/fctp->lumped_mass_matrix[i];            n_substeps=1.;
182              printf("%d : %e %e %e\n",i,u[i],fctp->u[i],rtmp);         }
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 critical
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                  /*
295  for (i=0;i<21;++i) {           * save last u
296  for (j=0;j<21;++j) printf("%d->%e ",i*21+j,fctp->u[i*21+j]);           *
297  printf("\n");           */
298            if (Paso_noError()) {
299                 #pragma omp parallel for schedule(static) private(i)
300                 for (i=0;i<n_rows;++i) {
301                     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  }  }
   
         t+=dt2;  
         n++;  
       }  
   }  
   if (Paso_noError()) {  
      #pragma omp parallel for schedule(static) private(i)  
      for (i=0;i<n_rows;++i) {  
        u[i]=fctp->u[i];  
      }  
   }  
 }        

Legend:
Removed from v.1388  
changed lines
  Added in v.1401

  ViewVC Help
Powered by ViewVC 1.1.26