proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
shockCapturing.c
Go to the documentation of this file.
1#include <math.h>
2#include "shockCapturing.h"
3
8void calculateNumericalDiffusionResGrad(int nElements_global,
9 int nQuadraturePoints_element,
10 int nSpace,
11 double shockCapturingDiffusion,
12 double* elementDiameter,
13 double* strong_residual,
14 double* grad_u,
15 double* numDiff)
16{
17 int eN,k,I;
18 double h,
19 num,
20 den,
21 n_grad_u[nQuadraturePoints_element],
22 n_R_eN,
23 maxNumDiff;
24 for(eN=0;eN<nElements_global;eN++)
25 {
26 h = elementDiameter[eN];
27 n_R_eN=0.0;
28 for (k=0;k<nQuadraturePoints_element;k++)
29 {
30 n_grad_u[k] = 0.0;
31 for (I=0;I<nSpace;I++)
32 {
33 n_grad_u[k] += grad_u[eN*nQuadraturePoints_element*nSpace+k*nSpace+I]*grad_u[eN*nQuadraturePoints_element*nSpace+k*nSpace+I];
34 }
35 n_grad_u[k] = sqrt(n_grad_u[k]);
36 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual[eN*nQuadraturePoints_element+k]);
37 den = n_grad_u[k];
38 if (den > num*1.0e-8)
39 numDiff[eN*nQuadraturePoints_element+k] = num/den;
40 else
41 numDiff[eN*nQuadraturePoints_element+k]=0.0;
42 n_R_eN=fmax(n_R_eN,fabs(strong_residual[eN*nQuadraturePoints_element+k]));
43 }
44 maxNumDiff=0.0;
45 for (k=0;k<nQuadraturePoints_element;k++)
46 {
47 num = shockCapturingDiffusion*0.5*h*n_R_eN;
48 den = n_grad_u[k];
49 if (den > num*1.0e-8)
50 numDiff[eN*nQuadraturePoints_element+k] = num/den;
51 else
52 numDiff[eN*nQuadraturePoints_element+k]=0.0;
53 maxNumDiff = fmax(maxNumDiff,numDiff[eN*nQuadraturePoints_element+k]);
54 }
55 for (k=0;k<nQuadraturePoints_element;k++)
56 numDiff[eN*nQuadraturePoints_element+k] = maxNumDiff;
57 }
58}
59
60void calculateNumericalDiffusionResGradQuad(int nElements_global,
61 int nQuadraturePoints_element,
62 int nSpace,
63 double shockCapturingDiffusion,
64 double* elementDiameter,
65 double* strong_residual,
66 double* grad_u,
67 double* numDiff)
68{
69 int eN,k,I;
70 double h,
71 num,
72 den,
73 n_grad_u;
74 for(eN=0;eN<nElements_global;eN++)
75 {
76 h = elementDiameter[eN];
77 for (k=0;k<nQuadraturePoints_element;k++)
78 {
79 n_grad_u = 0.0;
80 for (I=0;I<nSpace;I++)
81 {
82 n_grad_u += grad_u[eN*nQuadraturePoints_element*nSpace+k*nSpace+I]*grad_u[eN*nQuadraturePoints_element*nSpace+k*nSpace+I];
83 }
84 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual[eN*nQuadraturePoints_element+k]);
85 den = sqrt(n_grad_u) + 1.0e-8;
86 numDiff[eN*nQuadraturePoints_element+k] = num/den;
87 }
88 }
89}
90
91void calculateNumericalDiffusionEikonal(int nElements_global,
92 int nQuadraturePoints_element,
93 double shockCapturingDiffusion,
94 double* elementDiameter,
95 double* strong_residual,
96 double* numDiff)
97{
98 int eN,k;
99 double h;
100 for(eN=0;eN<nElements_global;eN++)
101 {
102 h = elementDiameter[eN];
103 for (k=0;k<nQuadraturePoints_element;k++)
104 {
105 numDiff[eN*nQuadraturePoints_element+k] = shockCapturingDiffusion*0.5*h*fabs(strong_residual[eN*nQuadraturePoints_element+k]);
106 //printf("%12.5e\n",numDiff[eN*nQuadraturePoints_element+k]);
107 }
108 }
109}
110
114void calculateNumericalDiffusionHJ(int nElements_global,
115 int nQuadraturePoints_element,
116 char shockCapturing,
117 double shockCapturingDiffusion,
118 double* elementDiameter,
119 double* strong_residual,
120 double* mt,
121 double* H,
122 double* numDiff)
123{
124 int eN,k;
125 double h,
126 num,
127 den,
128 n_R_eN;
129 //maxNumDiff;
130 for(eN=0;eN<nElements_global;eN++)
131 {
132 h = elementDiameter[eN];
133 n_R_eN=0.0;
134 for (k=0;k<nQuadraturePoints_element;k++)
135 {
136 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual[eN*nQuadraturePoints_element+k]);
137 den = fabs(H[eN*nQuadraturePoints_element+k]) +
138 fabs(mt[eN*nQuadraturePoints_element+k]);
139 if (den > num*1.0e-8)
140 numDiff[eN*nQuadraturePoints_element+k] = num/den;
141 else
142 numDiff[eN*nQuadraturePoints_element+k]=0.0;
143 n_R_eN=fmax(n_R_eN,fabs(strong_residual[eN*nQuadraturePoints_element+k]));
144 }
145/* maxNumDiff=0.0; */
146/* for (k=0;k<nQuadraturePoints_element;k++) */
147/* { */
148/* num = shockCapturingDiffusion*0.5*h*n_R_eN; */
149/* den = fabs(H[eN*nQuadraturePoints_element+k])+ */
150/* fabs(mt[eN*nQuadraturePoints_element+k]); */
151/* if (den > num*1.0e-8) */
152/* numDiff[eN*nQuadraturePoints_element+k] = num/den; */
153/* else */
154/* numDiff[eN*nQuadraturePoints_element+k]=0.0; */
155/* maxNumDiff = fmax(maxNumDiff,numDiff[eN*nQuadraturePoints_element+k]); */
156/* } */
157/* for (k=0;k<nQuadraturePoints_element;k++) */
158/* numDiff[eN*nQuadraturePoints_element+k] = maxNumDiff; */
159 }
160}
161
166void calculateNumericalDiffusionHJV2(int nElements_global,
167 int nQuadraturePoints_element,
168 char shockCapturing,
169 double shockCapturingDiffusion,
170 double* elementDiameter,
171 double* strong_residual,
172 double* mt,
173 double* H,
174 double* numDiff)
175{
176 int eN,k;
177 double h,num,den=0.0,enorm_h;
178 for(eN=0;eN<nElements_global;eN++)
179 {
180 h = elementDiameter[eN];
181 for (k=0;k<nQuadraturePoints_element;k++)
182 {
183 num = shockCapturingDiffusion*h*fabs(strong_residual[eN*nQuadraturePoints_element+
184 k]);
185 enorm_h = fabs(H[eN*nQuadraturePoints_element + k]);
186 if(shockCapturing == '1')
187 {
188 den = (fabs(mt[eN*nQuadraturePoints_element+
189 k]) +
190 enorm_h)+num*1.0e-8;
191 }
192 else if(shockCapturing == '2')
193 {
194 den = sqrt(mt[eN*nQuadraturePoints_element+
195 k]
196 *
197 mt[eN*nQuadraturePoints_element+
198 k]
199 +
200 enorm_h*enorm_h) + num*1.0e-8;
201 }
202 numDiff[eN*nQuadraturePoints_element+k] = num/den;
203 }
204 }
205}
206
207void calculateNumericalDiffusion_A_1(int nElements_global,
208 int nQuadraturePoints_element,
209 int nSpace,
210 double shockCapturingFactor,
211 double* elementDiameter,
212 double* strong_residual,
213 double* mt,
214 double* df,
215 double* numDiff)
216{
217 int eN,k,I;
218 double h,num,den,enorm_df;
219 for(eN=0;eN<nElements_global;eN++)
220 {
221 h = elementDiameter[eN];
222 for (k=0;k<nQuadraturePoints_element;k++)
223 {
224 num = shockCapturingFactor*h*fabs(strong_residual[eN*nQuadraturePoints_element+
225 k]);
226 enorm_df = 0.0;
227 for(I=0;I<nSpace;I++)
228 {
229 enorm_df
230 +=
231 df[eN*nQuadraturePoints_element*nSpace+
232 k*nSpace+
233 I]
234 *
235 df[eN*nQuadraturePoints_element*nSpace+
236 k*nSpace+
237 I];
238 }
239 enorm_df = sqrt(enorm_df);
240 den = (fabs(mt[eN*nQuadraturePoints_element+
241 k]) +
242 enorm_df);
243 numDiff[eN*nQuadraturePoints_element+k] = num/den;
244 }
245 }
246}
247/*mwf try Juanes' formula ...
248 D_{sc,g} = C_{sc}h|R|/|u/h| */
249void calculateNumericalDiffusionResGradJuanes(int nElements_global,
250 int nQuadraturePoints_element,
251 int nSpace,
252 double shockCapturingDiffusion,
253 double uSC,
254 double* elementDiameter,
255 double* strong_residual,
256 double* grad_u,
257 double* numDiff)
258{
259 int eN,k;
260 double h,
261 n_R_eN,
262 num,
263 den,
264 maxNumDiff;
265 for(eN=0;eN<nElements_global;eN++)
266 {
267 h = elementDiameter[eN];
268 n_R_eN=0.0;
269 for (k=0;k<nQuadraturePoints_element;k++)
270 {
271 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual[eN*nQuadraturePoints_element+k]);
272 den = fabs(uSC)/h;
273 if (den > num*1.0e-8)
274 numDiff[eN*nQuadraturePoints_element+k] = num/den;
275 else
276 numDiff[eN*nQuadraturePoints_element+k]=0.0;
277 n_R_eN=fmax(n_R_eN,fabs(strong_residual[eN*nQuadraturePoints_element+k]));
278 }
279 maxNumDiff=0.0;
280 for (k=0;k<nQuadraturePoints_element;k++)
281 {
282 num = shockCapturingDiffusion*0.5*h*n_R_eN;
283 den = fabs(uSC)/h;
284 if (den > num*1.0e-8)
285 numDiff[eN*nQuadraturePoints_element+k] = num/den;
286 else
287 numDiff[eN*nQuadraturePoints_element+k]=0.0;
288 maxNumDiff = fmax(maxNumDiff,numDiff[eN*nQuadraturePoints_element+k]);
289 }
290 for (k=0;k<nQuadraturePoints_element;k++)
291 numDiff[eN*nQuadraturePoints_element+k] = maxNumDiff;
292 }
293}
294
295/****
296 supposed to be nu_i = nu_c h^{2-beta}|R_i|
297 ***/
298
299void calculateNumericalDiffusionJaffre(int nElements_global,
300 int nQuadraturePoints_element,
301 int nSpace,
302 double shockCapturingDiffusion,
303 double beta,
304 double* elementDiameter,
305 double* strong_residual,
306 double* grad_u,
307 double* numDiff)
308{
309 int eN,k,I;
310 double h,
311 num,
312 hfactor,
313 n_grad_u[nQuadraturePoints_element],
314 n_R_eN,
315 maxNumDiff;
316
317 for(eN=0;eN<nElements_global;eN++)
318 {
319 h = elementDiameter[eN];
320 hfactor = pow(h,2.0-beta);
321 for (k=0;k<nQuadraturePoints_element;k++)
322 {
323 num = shockCapturingDiffusion*hfactor*fabs(strong_residual[eN*nQuadraturePoints_element+k]);
324 numDiff[eN*nQuadraturePoints_element+k] = num;
325 }
326
327 maxNumDiff=0.0;
328 for (k=0;k<nQuadraturePoints_element;k++)
329 {
330 maxNumDiff = fmax(maxNumDiff,numDiff[eN*nQuadraturePoints_element+k]);
331 }
332 for (k=0;k<nQuadraturePoints_element;k++)
333 numDiff[eN*nQuadraturePoints_element+k] = maxNumDiff;
334 }
335}
336
337
H
Double H
Definition Headers.h:65
num
Int num
Definition Headers.h:32
df
double df(double C, double b, double a, int q, int r)
Definition analyticalSolutions.c:2209
calculateNumericalDiffusionJaffre
void calculateNumericalDiffusionJaffre(int nElements_global, int nQuadraturePoints_element, int nSpace, double shockCapturingDiffusion, double beta, double *elementDiameter, double *strong_residual, double *grad_u, double *numDiff)
Definition shockCapturing.c:299
calculateNumericalDiffusionEikonal
void calculateNumericalDiffusionEikonal(int nElements_global, int nQuadraturePoints_element, double shockCapturingDiffusion, double *elementDiameter, double *strong_residual, double *numDiff)
Definition shockCapturing.c:91
calculateNumericalDiffusionResGradQuad
void calculateNumericalDiffusionResGradQuad(int nElements_global, int nQuadraturePoints_element, int nSpace, double shockCapturingDiffusion, double *elementDiameter, double *strong_residual, double *grad_u, double *numDiff)
Definition shockCapturing.c:60
calculateNumericalDiffusionHJV2
void calculateNumericalDiffusionHJV2(int nElements_global, int nQuadraturePoints_element, char shockCapturing, double shockCapturingDiffusion, double *elementDiameter, double *strong_residual, double *mt, double *H, double *numDiff)
Calculate the shock capturing diffusion for a Hamilton-Jacobi equation at the quadrature points mwf t...
Definition shockCapturing.c:166
calculateNumericalDiffusionResGradJuanes
void calculateNumericalDiffusionResGradJuanes(int nElements_global, int nQuadraturePoints_element, int nSpace, double shockCapturingDiffusion, double uSC, double *elementDiameter, double *strong_residual, double *grad_u, double *numDiff)
Definition shockCapturing.c:249
calculateNumericalDiffusionResGrad
void calculateNumericalDiffusionResGrad(int nElements_global, int nQuadraturePoints_element, int nSpace, double shockCapturingDiffusion, double *elementDiameter, double *strong_residual, double *grad_u, double *numDiff)
Definition shockCapturing.c:8
calculateNumericalDiffusionHJ
void calculateNumericalDiffusionHJ(int nElements_global, int nQuadraturePoints_element, char shockCapturing, double shockCapturingDiffusion, double *elementDiameter, double *strong_residual, double *mt, double *H, double *numDiff)
Calculate the shock capturing diffusion for a Hamilton-Jacobi equation at the quadrature points.
Definition shockCapturing.c:114
calculateNumericalDiffusion_A_1
void calculateNumericalDiffusion_A_1(int nElements_global, int nQuadraturePoints_element, int nSpace, double shockCapturingFactor, double *elementDiameter, double *strong_residual, double *mt, double *df, double *numDiff)
Definition shockCapturing.c:207
shockCapturing.h
C implementations of shock capturing diffusion calculations.
Generated on for proteus by doxygen 1.18.0