/[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 1369 by gross, Mon Dec 17 03:42:03 2007 UTC revision 1370 by gross, Wed Jan 2 09:21:43 2008 UTC
# Line 59  void Paso_FCTransportProblem_addAdvectiv Line 59  void Paso_FCTransportProblem_addAdvectiv
59             #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij)  schedule(static)             #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij)  schedule(static)
60             for (i = 0; i < n; ++i) {             for (i = 0; i < n; ++i) {
61                 if (fc->colorOf[i]==color) {                 if (fc->colorOf[i]==color) {
62                      fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]];
63    
64                    for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {                    for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
65                       j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];                       j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
66                       if (i<j) {                       if (j<i) {
67                          /* find entry a[j,i] */                          /* find entry a[j,i] */
68                          for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {                          for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
69                              if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {                              if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
70                                  d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],fc->flux_matrix->mainBlock->val[iptr_ji]);                                  d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],fc->flux_matrix->mainBlock->val[iptr_ji]);
71    printf("%d %d -> %e\n",i,j,d_ij);
72                                  fc->transport_matrix->mainBlock->val[iptr_ij]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;                                  fc->transport_matrix->mainBlock->val[iptr_ij]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;
73                                  fc->transport_matrix->mainBlock->val[iptr_ji]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ji]+d_ij;                                  fc->transport_matrix->mainBlock->val[iptr_ji]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ji]+d_ij;
74    printf("%d %d -> %e -> %e %e \n",i,j,d_ij,fc->transport_matrix->mainBlock->val[iptr_ij],fc->transport_matrix->mainBlock->val[iptr_ji]);
75                                  fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij;                                  fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij;
76                                  fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij;                                  fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij;
77                                  break;                                  break;
# Line 84  void Paso_FCTransportProblem_addAdvectiv Line 88  void Paso_FCTransportProblem_addAdvectiv
88    }    }
89    
90  }  }
91    
92    void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {
93    
94      double *remote_u=NULL;
95      
96      if (fc==NULL) return;
97      Paso_SystemMatrix_startCollect(fc->transport_matrix,u);
98      /* process main block */
99      Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa);
100      /* finish exchange */
101      remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix);
102      /* process couple block */
103      Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa);
104    
105      Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa);
106    }
107  /**************************************************************/  /**************************************************************/
108    
109  /* adds antidiffusion to fa    /* adds antidiffusion to fa  
# Line 102  void Paso_FCTransportProblem_addAdvectiv Line 122  void Paso_FCTransportProblem_addAdvectiv
122    
123  */  */
124    
125  void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {  void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {
126    
127    register double u_i,P_p,P_n,Q_p,Q_n,r_p,r_n, a_ij, d, u_j, r_ij, f_ij, a_ji;    register double u_i,P_p,P_n,Q_p,Q_n,r_p,r_n, a_ij, d, u_j, r_ij, f_ij, a_ji;
   double *u_remote=NULL;  
128    index_t color, iptr_ij,j,iptr_ji, i;    index_t color, iptr_ij,j,iptr_ji, i;
129    dim_t n;    dim_t n;
130    
# Line 113  void Paso_FCTransportProblem_setAntiDiff Line 132  void Paso_FCTransportProblem_setAntiDiff
132    if (fc==NULL) return;    if (fc==NULL) return;
133    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
134    /* exchange */    /* exchange */
135    Paso_SystemMatrix_allocBuffer(fc->flux_matrix);  
   Paso_SystemMatrix_startCollect(fc->flux_matrix,u);  
   u_remote=Paso_SystemMatrix_finishCollect(fc->flux_matrix);  
136    
137    #pragma omp parallel private(color)    #pragma omp parallel private(color)
138    {    {

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

  ViewVC Help
Powered by ViewVC 1.1.26