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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1361 - (hide annotations)
Fri Dec 14 09:26:51 2007 UTC (12 years, 10 months ago) by gross
Original Path: trunk/paso/src/Solver_FluxControl.c
File MIME type: text/plain
File size: 10444 byte(s)
first steps towards a flux controlled transport solver.
1 gross 1361 /* $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(2*ABS(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     /* free all memory used by FluxControl */
40    
41     void Paso_Solver_FluxControl_free(Paso_Solver_FluxControl* in) {
42     if (in!=NULL) {
43     Paso_SystemMatrix_freeBuffer(in->matrix);
44     Paso_SystemMatrix_free(in->matrix);
45     MEMFREE(in->colorOf);
46     MEMFREE(in->main_iptr);
47     MEMFREE(in);
48     }
49     }
50    
51     /**************************************************************/
52    
53     /* constructs a flux control mechanism */
54    
55     Paso_Solver_FluxControl* Paso_SolverFCT_getFluxControl(Paso_SystemMatrix * A) {
56    
57     Paso_Solver_FluxControl* out=NULL;
58     dim_t n,i;
59     index_t iptr,iptr_main,k;
60    
61     if (A==NULL) return out;
62     n=Paso_SystemMatrix_getTotalNumRows(A);
63     if (A->block_size!=1) {
64     Paso_setError(TYPE_ERROR,"Paso_SolverFCT_getFluxControl: block size > 1 is not supported.");
65     return NULL;
66     }
67     out=MEMALLOC(1,Paso_Solver_FluxControl);
68     if (Paso_checkPtr(out)) return NULL;
69    
70     out->matrix=Paso_SystemMatrix_reference(A);
71     out->colorOf=NULL;
72     out->main_iptr=NULL;
73    
74    
75     /* allocations: */
76     out->colorOf=MEMALLOC(n,index_t);
77     out->main_iptr=MEMALLOC(n,index_t);
78     if ( ! (Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) ) ) {
79     printf("Paso_SolverFCT_getFluxControl: Revise coloring!!\n");
80     Paso_Pattern_color(A->mainBlock->pattern,&(out->num_colors),out->colorOf);
81     Paso_SystemMatrix_allocBuffer(A);
82    
83     #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)
84     for (i = 0; i < n; ++i) {
85     for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; ++iptr) {
86     iptr_main=A->mainBlock->pattern->ptr[0]-1;
87     for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; iptr++) {
88     if (A->mainBlock->pattern->index[iptr]==i) {
89     iptr_main=iptr;
90     break;
91     }
92     }
93     out->main_iptr[i]=iptr_main;
94     if (iptr_main==A->mainBlock->pattern->ptr[0]-1)
95     Paso_setError(VALUE_ERROR, "Paso_SolverFCT_getFluxControl: no main diagonal");
96     }
97     }
98    
99     }
100     if (Paso_noError()) {
101     return out;
102     } else {
103     Paso_Solver_FluxControl_free(out);
104     return NULL;
105     }
106     }
107    
108     /**************************************************************/
109    
110     /* adds A plus stabelising diffusion into the matrix B */
111    
112     /* d_ij=alpha*max(0,-a[i,j],-a[j,i]) */
113     /* b[i,j]+=alpha*(a[i,j]+d_ij) */
114     /* b[j,i]+=alpha*(a[j,i]+d_ij) */
115     /* b[i,i]-=alpha*d_ij */
116     /* b[j,j]-=alpha*d_ij */
117    
118     void Paso_Solver_FluxControl_addDiffusion(Paso_Solver_FluxControl * fc, double alpha, Paso_SystemMatrix * B) {
119     dim_t n,i;
120     index_t color, iptr_ij,j,iptr_ji;
121     register double d_ij;
122    
123     if (fc==NULL) return;
124     n=Paso_SystemMatrix_getTotalNumRows(fc->matrix);
125     /* TODO test - same pattern + block size */
126    
127     #pragma omp parallel private(color)
128     {
129     /* process main block */
130     for (color=0;color<fc->num_colors;++color) {
131     #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij) schedule(static)
132     for (i = 0; i < n; ++i) {
133     if (fc->colorOf[i]==color) {
134     for (iptr_ij=fc->matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
135     j=fc->matrix->mainBlock->pattern->index[iptr_ij];
136     if (i<j) {
137     /* find entry a[j,i] */
138     for (iptr_ji=fc->matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {
139     if (fc->matrix->mainBlock->pattern->index[iptr_ji]==i) {
140     d_ij=(-alpha)*MIN3(0.,fc->matrix->mainBlock->val[iptr_ij],fc->matrix->mainBlock->val[iptr_ji]);
141     B->mainBlock->val[iptr_ij]+=alpha*fc->matrix->mainBlock->val[iptr_ij]+d_ij;
142     B->mainBlock->val[iptr_ji]+=alpha*fc->matrix->mainBlock->val[iptr_ji]+d_ij;
143     B->mainBlock->val[fc->main_iptr[i]]-=d_ij;
144     B->mainBlock->val[fc->main_iptr[j]]-=d_ij;
145     break;
146     }
147     }
148     }
149    
150     }
151     /* TODO process couple block */
152     }
153     }
154     #pragma omp barrier
155     }
156     }
157    
158     }
159     /**************************************************************/
160    
161     /* adds antidiffusion to fa
162    
163     P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i])
164     P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i])
165     Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i])
166     Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i])
167     d_ij=max(0,-a[i,j],-a[j,i])
168     l_ji=max(0,a[j,i],a[j,i]-a[i,j])
169     if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij)
170     r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i]
171     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])
172     fa[i]+=f_ij
173     fa[j]-=f_ij
174    
175     */
176    
177     void Paso_Solver_FluxControl_setAntiDiffusiveFlux(Paso_Solver_FluxControl * fc, double * u, double* fa) {
178    
179     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;
180     double *u_remote=NULL;
181     index_t color, iptr_ij,j,iptr_ji, i;
182     dim_t n;
183    
184    
185     if (fc==NULL) return;
186     n=Paso_SystemMatrix_getTotalNumRows(fc->matrix);
187     /* exchange */
188     Paso_SystemMatrix_startCollect(fc->matrix,u);
189     u_remote=Paso_SystemMatrix_finishCollect(fc->matrix);
190    
191     #pragma omp parallel private(color)
192     {
193     for (color=0;color<fc->num_colors;++color) {
194     #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)
195     for (i = 0; i < n; ++i) {
196     if (fc->colorOf[i]==color) {
197     u_i=u[i];
198     /* gather the smoothness sensor */
199     P_p=0.;
200     P_n=0.;
201     Q_p=0.;
202     Q_n=0.;
203     #pragma ivdep
204     for (iptr_ij=(fc->matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {
205     a_ij=fc->matrix->mainBlock->val[iptr_ij];
206     j=fc->matrix->mainBlock->pattern->index[iptr_ij];
207     d=u[j]-u_i;
208     if (a_ij<0.) {
209     if (d<0.) {
210     P_p+=a_ij*d;
211     } else {
212     P_n+=a_ij*d;
213     }
214     } else {
215     if (d>0.) {
216     Q_p+=a_ij*d;
217     } else {
218     Q_n+=a_ij*d;
219     }
220     }
221     }
222     #pragma ivdep
223     for (iptr_ij=(fc->matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {
224     a_ij=fc->matrix->coupleBlock->val[iptr_ij];
225     j=fc->matrix->coupleBlock->pattern->index[iptr_ij];
226     d=u_remote[j]-u_i;
227     if (a_ij<0.) {
228     if (d<0.) {
229     P_p+=a_ij*d;
230     } else {
231     P_n+=a_ij*d;
232     }
233     } else {
234     if (d>0.) {
235     Q_p+=a_ij*d;
236     } else {
237     Q_n+=a_ij*d;
238     }
239     }
240     }
241     /* set the smoothness indicators */
242     r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
243     r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.;
244     /* anti diffusive flux from main block */
245     for (iptr_ij=fc->matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
246     a_ij=fc->matrix->mainBlock->val[iptr_ij];
247     j=fc->matrix->mainBlock->pattern->index[iptr_ij];
248     if (a_ij < 0 && i!=j) {
249     /* find entry a[j,i] */
250     for (iptr_ji=fc->matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {
251     if (fc->matrix->mainBlock->pattern->index[iptr_ji]==i) {
252     a_ji=fc->matrix->mainBlock->val[iptr_ji];
253     if (a_ji > a_ij || (a_ji == a_ij && j<i) ) {
254     u_j=u[j];
255     r_ij = u_i>u_j ? r_p : r_n;
256     f_ij =MIN(r_ij*a_ij,a_ji-a_ij)*(u_i-u_j);
257     fa[i]+=f_ij;
258     fa[j]-=f_ij;
259     break;
260     }
261     }
262     }
263     }
264     }
265     /* anti diffusive flux from couple block */
266    
267     /* TODO */
268     }
269     }
270     #pragma omp barrier
271     }
272     }
273     }

  ViewVC Help
Powered by ViewVC 1.1.26