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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1362 - (show annotations)
Mon Dec 17 02:28:16 2007 UTC (11 years, 4 months ago) by gross
File MIME type: text/plain
File size: 8293 byte(s)
and more on FCT solver
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(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 /* adds A plus stabelising diffusion into the matrix B */
40
41 /* d_ij=alpha*max(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]-=alpha*d_ij */
45 /* b[j,j]-=alpha*d_ij */
46
47 void Paso_FCTransportProblem_addDiffusion(Paso_FCTransportProblem * fc, double alpha, Paso_SystemMatrix * B) {
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 /* TODO test - same pattern + block size */
55
56 #pragma omp parallel private(color)
57 {
58 /* process main block */
59 for (color=0;color<fc->num_colors;++color) {
60 #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij) schedule(static)
61 for (i = 0; i < n; ++i) {
62 if (fc->colorOf[i]==color) {
63 for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
64 j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
65 if (i<j) {
66 /* find entry a[j,i] */
67 for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {
68 if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
69 d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],fc->flux_matrix->mainBlock->val[iptr_ji]);
70 B->mainBlock->val[iptr_ij]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;
71 B->mainBlock->val[iptr_ji]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ji]+d_ij;
72 B->mainBlock->val[fc->main_iptr[i]]-=d_ij;
73 B->mainBlock->val[fc->main_iptr[j]]-=d_ij;
74 break;
75 }
76 }
77 }
78
79 }
80 /* TODO process couple block */
81 }
82 }
83 #pragma omp barrier
84 }
85 }
86
87 }
88 /**************************************************************/
89
90 /* adds antidiffusion to fa
91
92 P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i])
93 P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i])
94 Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i])
95 Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i])
96 d_ij=max(0,-a[i,j],-a[j,i])
97 l_ji=max(0,a[j,i],a[j,i]-a[i,j])
98 if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij)
99 r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i]
100 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])
101 fa[i]+=f_ij
102 fa[j]-=f_ij
103
104 */
105
106 void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {
107
108 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;
109 double *u_remote=NULL;
110 index_t color, iptr_ij,j,iptr_ji, i;
111 dim_t n;
112
113
114 if (fc==NULL) return;
115 n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
116 /* exchange */
117 Paso_SystemMatrix_allocBuffer(fc->flux_matrix);
118 Paso_SystemMatrix_startCollect(fc->flux_matrix,u);
119 u_remote=Paso_SystemMatrix_finishCollect(fc->flux_matrix);
120
121 #pragma omp parallel private(color)
122 {
123 for (color=0;color<fc->num_colors;++color) {
124 #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)
125 for (i = 0; i < n; ++i) {
126 if (fc->colorOf[i]==color) {
127 u_i=u[i];
128 /* gather the smoothness sensor */
129 P_p=0.;
130 P_n=0.;
131 Q_p=0.;
132 Q_n=0.;
133 #pragma ivdep
134 for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {
135 a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
136 j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
137 d=u[j]-u_i;
138 if (a_ij<0.) {
139 if (d<0.) {
140 P_p+=a_ij*d;
141 } else {
142 P_n+=a_ij*d;
143 }
144 } else {
145 if (d>0.) {
146 Q_p+=a_ij*d;
147 } else {
148 Q_n+=a_ij*d;
149 }
150 }
151 }
152 #pragma ivdep
153 for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {
154 a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];
155 j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];
156 d=u_remote[j]-u_i;
157 if (a_ij<0.) {
158 if (d<0.) {
159 P_p+=a_ij*d;
160 } else {
161 P_n+=a_ij*d;
162 }
163 } else {
164 if (d>0.) {
165 Q_p+=a_ij*d;
166 } else {
167 Q_n+=a_ij*d;
168 }
169 }
170 }
171 /* set the smoothness indicators */
172 r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
173 r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.;
174 /* anti diffusive flux from main block */
175 for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
176 a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
177 j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
178 if (a_ij < 0 && i!=j) {
179 /* find entry a[j,i] */
180 for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {
181 if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
182 a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];
183 if (a_ji > a_ij || (a_ji == a_ij && j<i) ) {
184 u_j=u[j];
185 r_ij = u_i>u_j ? r_p : r_n;
186 f_ij =MIN(r_ij*a_ij,a_ji-a_ij)*(u_i-u_j);
187 fa[i]+=f_ij;
188 fa[j]-=f_ij;
189 break;
190 }
191 }
192 }
193 }
194 }
195 /* anti diffusive flux from couple block */
196
197 /* TODO */
198 }
199 }
200 #pragma omp barrier
201 }
202 }
203 }

  ViewVC Help
Powered by ViewVC 1.1.26