/[escript]/trunk/paso/src/SolverFCT_FluxControl.c
ViewVC logotype

Diff of /trunk/paso/src/SolverFCT_FluxControl.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1408 by gross, Mon Feb 4 07:19:50 2008 UTC revision 1410 by gross, Thu Feb 7 04:24:00 2008 UTC
# Line 163  void Paso_SolverFCT_setQs(const double* Line 163  void Paso_SolverFCT_setQs(const double*
163           j=pattern->mainPattern->index[iptr_ij];           j=pattern->mainPattern->index[iptr_ij];
164           u_j=u[j];           u_j=u[j];
165           u_min_i=MIN(u_min_i,u_j);           u_min_i=MIN(u_min_i,u_j);
166           u_max_i=MIN(u_max_i,u_j);           u_max_i=MAX(u_max_i,u_j);
167       }       }
168       #pragma ivdep       #pragma ivdep
169       for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {       for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
170            j=pattern->couplePattern->index[iptr_ij];            j=pattern->couplePattern->index[iptr_ij];
171            u_j=remote_u[j];            u_j=remote_u[j];
172            u_min_i=MIN(u_min_i,u_j);            u_min_i=MIN(u_min_i,u_j);
173            u_max_i=MIN(u_max_i,u_j);            u_max_i=MAX(u_max_i,u_j);
174        }        }
175        QN[i]=u_min_i-u_i;        QN[i]=u_min_i-u_i;
176        QP[i]=u_max_i-u_i;        QP[i]=u_max_i-u_i;
# Line 179  void Paso_SolverFCT_setQs(const double* Line 179  void Paso_SolverFCT_setQs(const double*
179    
180  /*  /*
181   *   *
182     *  f_{ij} = (m_{ij} - dt (1-theta) d_{ij}) (u_last[j]-u_last[i]) - (m_{ij} + dt theta d_{ij}) (u[j]-u[i])
183     *        
184     * m=fc->mass matrix
185     * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)
186     */
187    void Paso_FCTransportProblem_setAntiDiffusionFlux(const double dt, const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix, const double* u, const double* u_last)
188    {
189      dim_t n, j, i;
190      double *remote_u=NULL, *remote_u_last=NULL;
191      index_t iptr_ij;
192      const double f1=- dt * (1.-fc->theta);
193      const double f2=  dt * fc->theta;
194      register double m_ij, d_ij, u_i, u_last_i, d_u_last, d_u;
195      Paso_SystemMatrixPattern *pattern;
196      Paso_SystemMatrix_startCollect(fc->iteration_matrix,u);
197      Paso_SystemMatrix_startCollect(fc->iteration_matrix,u_last);
198      n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
199      pattern=fc->iteration_matrix->pattern;
200      remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);
201      remote_u_last=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);
202      if ( (ABS(f1) >0 ) ) {
203         if ( (ABS(f2) >0 ) ) {
204            #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
205            for (i = 0; i < n; ++i) {
206               u_i=u[i];
207               u_last_i=u_last[i];
208               #pragma ivdep
209               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
210                  j=pattern->mainPattern->index[iptr_ij];
211                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
212                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
213                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
214                  d_u=u[j]-u_i;
215                  d_u_last=u_last[j]-u_last_i;
216                  flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
217               }
218               #pragma ivdep
219               for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
220                  j=pattern->couplePattern->index[iptr_ij];
221                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
222                  d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
223                                                 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
224                  d_u=remote_u[j]-u_i;
225                  d_u_last=remote_u_last[j]-u_last_i;
226                  flux_matrix->coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
227               }
228            }
229         } else {
230            #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
231            for (i = 0; i < n; ++i) {
232               u_i=u[i];
233               u_last_i=u_last[i];
234               #pragma ivdep
235               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
236                  j=pattern->mainPattern->index[iptr_ij];
237                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
238                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
239                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
240                  d_u=u[j]-u_i;
241                  d_u_last=u_last[j]-u_last_i;
242                  flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
243               }
244               #pragma ivdep
245               for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
246                  j=pattern->couplePattern->index[iptr_ij];
247                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
248                  d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
249                                                 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
250                  d_u=remote_u[j]-u_i;
251                  d_u_last=remote_u_last[j]-u_last_i;
252                  flux_matrix->coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
253               }
254            }
255         }
256      } else {
257         if ( (ABS(f2) >0 ) ) {
258            #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
259            for (i = 0; i < n; ++i) {
260               u_i=u[i];
261               u_last_i=u_last[i];
262               #pragma ivdep
263               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
264                  j=pattern->mainPattern->index[iptr_ij];
265                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
266                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
267                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
268                  d_u=u[j]-u_i;
269                  d_u_last=u_last[j]-u_last_i;
270                  flux_matrix->mainBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
271               }
272               #pragma ivdep
273               for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
274                  j=pattern->couplePattern->index[iptr_ij];
275                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
276                  d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
277                                                 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
278                  d_u=remote_u[j]-u_i;
279                  d_u_last=remote_u_last[j]-u_last_i;
280                  flux_matrix->coupleBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
281               }
282            }
283         } else {
284            #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
285            for (i = 0; i < n; ++i) {
286               u_i=u[i];
287               u_last_i=u_last[i];
288               #pragma ivdep
289               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
290                  j=pattern->mainPattern->index[iptr_ij];
291                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
292                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
293                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
294                  d_u=u[j]-u_i;
295                  d_u_last=u_last[j]-u_last_i;
296                  flux_matrix->mainBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
297               }
298               #pragma ivdep
299               for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
300                  j=pattern->couplePattern->index[iptr_ij];
301                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
302                  d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
303                                                 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
304                  d_u=remote_u[j]-u_i;
305                  d_u_last=remote_u_last[j]-u_last_i;
306                  flux_matrix->coupleBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
307               }
308            }
309          }
310      }
311    }
312    /*
313     *
314   * f_{ij} + = (a*m_{ij} + b* d_{ij}) (u[j]-u[i])   * f_{ij} + = (a*m_{ij} + b* d_{ij}) (u[j]-u[i])
315   *   *
316   * m=fc->mass matrix   * m=fc->mass matrix
# Line 196  void Paso_FCTransportProblem_updateAntiD Line 328  void Paso_FCTransportProblem_updateAntiD
328    n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);    n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
329    pattern=fc->iteration_matrix->pattern;    pattern=fc->iteration_matrix->pattern;
330    remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);    remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);
   
331    if ( (ABS(a) >0 ) ) {    if ( (ABS(a) >0 ) ) {
332        if ( (ABS(b) >0 ) ) {        if ( (ABS(b) >0 ) ) {
333           #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij,k_ij,ml_ij)           #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij,k_ij,ml_ij)
# Line 316  void Paso_FCTransportProblem_setRs(const Line 447  void Paso_FCTransportProblem_setRs(const
447            j=pattern->mainPattern->index[iptr_ij];            j=pattern->mainPattern->index[iptr_ij];
448            if (i != j ) {            if (i != j ) {
449               f_ij=f->mainBlock->val[iptr_ij];               f_ij=f->mainBlock->val[iptr_ij];
450               if (f_ij <=0  ) {               if (f_ij <=0) {
451                  PN_i+=f_ij;                  PN_i+=f_ij;
452               } else {               } else {
453                  PP_i+=f_ij;                  PP_i+=f_ij;
# Line 382  void Paso_FCTransportProblem_addCorrecte Line 513  void Paso_FCTransportProblem_addCorrecte
513                f_i+=f_ij*MIN(RP_i,remote_RN[j]);                f_i+=f_ij*MIN(RP_i,remote_RN[j]);
514            }            }
515        }        }
516        f[i]=f_i;        f[i]+=f_i;
517    }    }
518    Paso_SystemMatrix_startCollect(flux_matrix,RP);    Paso_SystemMatrix_startCollect(flux_matrix,RP);
519    remote_RP=Paso_SystemMatrix_finishCollect(flux_matrix);    remote_RP=Paso_SystemMatrix_finishCollect(flux_matrix);

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

  ViewVC Help
Powered by ViewVC 1.1.26