paso/src/SolverFCT_FluxControl.c

/* $Id:$ */

/*******************************************************
 *
 *       Copyright 2007 by University of Queensland
 *
 *                http://esscc.uq.edu.au
 *        Primary Business: Queensland, Australia
 *  Licensed under the Open Software License version 3.0
 *     http://www.opensource.org/licenses/osl-3.0.php
 *
 *******************************************************/

/**************************************************************/

/* Paso: FluxControl                                          */

/**************************************************************/

/* Author: l.gross@uq.edu.au                                */

/**************************************************************/


#include "Paso.h"
#include "SolverFCT.h"
#include "PasoUtil.h"

#define FLUX_S(a,b) ((SIGN(a)+SIGN(b))/2.)
#define MINMOD(a,b) (FLUX_S(a,b)*MIN(ABS(a),ABS(b)))
#define SUPERBEE(a,b) (FLUX_S(a,b)*MAX(MIN(2*ABS(a),ABS(b)),MIN(ABS(a),2*ABS(b))))

#define FLUX_L(a,b) SUPERBEE(a,b)  /* alter for other flux limiter */

#define FLUX_LIMITER(a) FLUX_L(a,1)

/**************************************************************/

/* adds A plus stabelising diffusion into the matrix B        */

/* d_ij=alpha*max(0,-a[i,j],-a[j,i])  */
/* b[i,j]+=alpha*(a[i,j]+d_ij)  */
/* b[j,i]+=alpha*(a[j,i]+d_ij)  */
/* b[i,i]-=alpha*d_ij  */
/* b[j,j]-=alpha*d_ij  */

void Paso_FCTransportProblem_addAdvectivePart(Paso_FCTransportProblem * fc, double alpha) {
  dim_t n,i;
  index_t color, iptr_ij,j,iptr_ji;
  register double d_ij;

  if (fc==NULL) return;
  n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);

  #pragma omp parallel private(color) 
  {
       /* process main block */
       for (color=0;color<fc->num_colors;++color) {
           #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij)  schedule(static)
           for (i = 0; i < n; ++i) {
               if (fc->colorOf[i]==color) {
                  fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]];

                  for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
                     j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
                     if (j<i) {
                        /* find entry a[j,i] */
                        for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
                            if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
                                d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],fc->flux_matrix->mainBlock->val[iptr_ji]);
printf("%d %d -> %e\n",i,j,d_ij);
                                fc->transport_matrix->mainBlock->val[iptr_ij]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;
                                fc->transport_matrix->mainBlock->val[iptr_ji]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ji]+d_ij;
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]);
                                fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij;
                                fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij;
                                break;
                            }
                        }
                     }
                        
                  }
                  /* TODO process couple block */
               }
           }
           #pragma omp barrier
       }
  }

}

void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {

  double *remote_u=NULL;
  
  if (fc==NULL) return;
  Paso_SystemMatrix_startCollect(fc->transport_matrix,u);
  /* process main block */
  Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa);
  /* finish exchange */
  remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix);
  /* process couple block */
  Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa);

  Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa); 
}
/**************************************************************/

/* adds antidiffusion to fa  
 
   P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i])   
   P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i])  
   Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i])  
   Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i])   
   d_ij=max(0,-a[i,j],-a[j,i])
   l_ji=max(0,a[j,i],a[j,i]-a[i,j])
   if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij)
      r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i]
      f_ij =min(FLUX_LIMITER(r_ij)*d_ij,l_ji) (u[i]-u[j])=min(FLUX_LIMITER(r_ij)*a[i,j],a[j,i]-a[i,j]) (u[i]-u[j])
      fa[i]+=f_ij
      fa[j]-=f_ij

*/

void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {

  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;
  index_t color, iptr_ij,j,iptr_ji, i;
  dim_t n;


  if (fc==NULL) return;
  n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
  /* exchange */


  #pragma omp parallel private(color) 
  {
       for (color=0;color<fc->num_colors;++color) {
           #pragma omp for schedule(static) private(i, u_i,P_p,P_n,Q_p,Q_n,r_p,r_n,iptr_ij,a_ij,d,j,iptr_ji, u_j, r_ij, f_ij, a_ji)
           for (i = 0; i < n; ++i) {
              if (fc->colorOf[i]==color) {
                  u_i=u[i];
                  /* gather the smoothness sensor */
                  P_p=0.;
                  P_n=0.;
                  Q_p=0.;
                  Q_n=0.;
                  #pragma ivdep
              for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {
                      a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
                      j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
                      d=u[j]-u_i;
                      if (a_ij<0.) {
                         if (d<0.) {
                            P_p+=a_ij*d;
                         } else {
                            P_n+=a_ij*d;
                         }
                      } else {
                         if (d>0.) {
                            Q_p+=a_ij*d;
                         } else {
                            Q_n+=a_ij*d;
                         }
                      }
              }
                  #pragma ivdep
              for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {
                      a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];
                      j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];
                      d=u_remote[j]-u_i;
                      if (a_ij<0.) {
                         if (d<0.) {
                            P_p+=a_ij*d;
                         } else {
                            P_n+=a_ij*d;
                         }
                      } else {
                         if (d>0.) {
                            Q_p+=a_ij*d;
                         } else {
                            Q_n+=a_ij*d;
                         }
                      }
              }
                  /* set the smoothness indicators */
                  r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
                  r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.;
                  /* anti diffusive flux from main block */
                  for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
                     a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
                     j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
                     if (a_ij < 0 && i!=j) {
                        /* find entry a[j,i] */
                        for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {
                            if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
                                a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];
                                if  (a_ji > a_ij || (a_ji == a_ij && j<i) ) {
                                    u_j=u[j];
                                    r_ij = u_i>u_j ? r_p : r_n;
                                    f_ij =MIN(r_ij*a_ij,a_ji-a_ij)*(u_i-u_j);
                                    fa[i]+=f_ij;
                                    fa[j]-=f_ij;
                                    break;
                                }
                            }
                        }
                     }
                  }
                  /* anti diffusive flux from couple block */

                  /* TODO */
              }
           }
           #pragma omp barrier
       }
  }
}
1	/* $Id:$ */
2
3	/*******************************************************
4	*
5	* Copyright 2007 by University of Queensland
6	*
7	* http://esscc.uq.edu.au
8	* Primary Business: Queensland, Australia
9	* Licensed under the Open Software License version 3.0
10	* http://www.opensource.org/licenses/osl-3.0.php
11	*
12	*******************************************************/
13
14	/**************************************************************/
15
16	/* Paso: FluxControl */
17
18	/**************************************************************/
19
20	/* Author: l.gross@uq.edu.au */
21
22	/**************************************************************/
23
24
25	#include "Paso.h"
26	#include "SolverFCT.h"
27	#include "PasoUtil.h"
28
29	#define FLUX_S(a,b) ((SIGN(a)+SIGN(b))/2.)
30	#define MINMOD(a,b) (FLUX_S(a,b)*MIN(ABS(a),ABS(b)))
31	#define SUPERBEE(a,b) (FLUX_S(a,b)MAX(MIN(2ABS(a),ABS(b)),MIN(ABS(a),2*ABS(b))))
32
33	#define FLUX_L(a,b) SUPERBEE(a,b) /* alter for other flux limiter */
34
35	#define FLUX_LIMITER(a) FLUX_L(a,1)
36
37	/**************************************************************/
38
39	/* adds A plus stabelising diffusion into the matrix B */
40
41	/* d_ij=alphamax(0,-a[i,j],-a[j,i]) /
42	/* b[i,j]+=alpha(a[i,j]+d_ij) /
43	/* b[j,i]+=alpha(a[j,i]+d_ij) /
44	/* b[i,i]-=alphad_ij /
45	/* b[j,j]-=alphad_ij /
46
47	void Paso_FCTransportProblem_addAdvectivePart(Paso_FCTransportProblem * fc, double alpha) {
48	dim_t n,i;
49	index_t color, iptr_ij,j,iptr_ji;
50	register double d_ij;
51
52	if (fc==NULL) return;
53	n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
54
55	#pragma omp parallel private(color)
56	{
57	/* process main block */
58	for (color=0;color<fc->num_colors;++color) {
59	#pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij) schedule(static)
60	for (i = 0; i < n; ++i) {
61	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) {
65	j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
66	if (j<i) {
67	/* find entry a[j,i] */
68	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) {
70	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;
73	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;
76	fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij;
77	break;
78	}
79	}
80	}
81
82	}
83	/* TODO process couple block */
84	}
85	}
86	#pragma omp barrier
87	}
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
110
111	P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i])
112	P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i])
113	Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i])
114	Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i])
115	d_ij=max(0,-a[i,j],-a[j,i])
116	l_ji=max(0,a[j,i],a[j,i]-a[i,j])
117	if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij)
118	r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i]
119	f_ij =min(FLUX_LIMITER(r_ij)d_ij,l_ji) (u[i]-u[j])=min(FLUX_LIMITER(r_ij)a[i,j],a[j,i]-a[i,j]) (u[i]-u[j])
120	fa[i]+=f_ij
121	fa[j]-=f_ij
122
123	*/
124
125	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;
128	index_t color, iptr_ij,j,iptr_ji, i;
129	dim_t n;
130
131
132	if (fc==NULL) return;
133	n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
134	/* exchange */
135
136
137	#pragma omp parallel private(color)
138	{
139	for (color=0;color<fc->num_colors;++color) {
140	#pragma omp for schedule(static) private(i, u_i,P_p,P_n,Q_p,Q_n,r_p,r_n,iptr_ij,a_ij,d,j,iptr_ji, u_j, r_ij, f_ij, a_ji)
141	for (i = 0; i < n; ++i) {
142	if (fc->colorOf[i]==color) {
143	u_i=u[i];
144	/* gather the smoothness sensor */
145	P_p=0.;
146	P_n=0.;
147	Q_p=0.;
148	Q_n=0.;
149	#pragma ivdep
150	for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {
151	a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
152	j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
153	d=u[j]-u_i;
154	if (a_ij<0.) {
155	if (d<0.) {
156	P_p+=a_ij*d;
157	} else {
158	P_n+=a_ij*d;
159	}
160	} else {
161	if (d>0.) {
162	Q_p+=a_ij*d;
163	} else {
164	Q_n+=a_ij*d;
165	}
166	}
167	}
168	#pragma ivdep
169	for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {
170	a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];
171	j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];
172	d=u_remote[j]-u_i;
173	if (a_ij<0.) {
174	if (d<0.) {
175	P_p+=a_ij*d;
176	} else {
177	P_n+=a_ij*d;
178	}
179	} else {
180	if (d>0.) {
181	Q_p+=a_ij*d;
182	} else {
183	Q_n+=a_ij*d;
184	}
185	}
186	}
187	/* set the smoothness indicators */
188	r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
189	r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.;
190	/* anti diffusive flux from main block */
191	for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
192	a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
193	j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
194	if (a_ij < 0 && i!=j) {
195	/* find entry a[j,i] */
196	for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {
197	if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
198	a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];
199	if (a_ji > a_ij \|\| (a_ji == a_ij && j<i) ) {
200	u_j=u[j];
201	r_ij = u_i>u_j ? r_p : r_n;
202	f_ij =MIN(r_ija_ij,a_ji-a_ij)(u_i-u_j);
203	fa[i]+=f_ij;
204	fa[j]-=f_ij;
205	break;
206	}
207	}
208	}
209	}
210	}
211	/* anti diffusive flux from couple block */
212
213	/* TODO */
214	}
215	}
216	#pragma omp barrier
217	}
218	}
219	}