proteus  1.8.1
C/C++/Fortran libraries
shockCapturing.c
Go to the documentation of this file.
1 #include <math.h>
2 #include "shockCapturing.h"
3 
8 void 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 
60 void 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 
91 void 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 
114 void 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 
166 void 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 
207 void 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| */
249 void 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 
299 void 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 
df
double df(double C, double b, double a, int q, int r)
Definition: analyticalSolutions.c:2209
calculateNumericalDiffusionEikonal
void calculateNumericalDiffusionEikonal(int nElements_global, int nQuadraturePoints_element, double shockCapturingDiffusion, double *elementDiameter, double *strong_residual, double *numDiff)
Definition: shockCapturing.c:91
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
num
Int num
Definition: Headers.h:32
H
Double H
Definition: Headers.h:65
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
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
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
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
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
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
shockCapturing.h
C implementations of shock capturing diffusion calculations.
Generated on Fri Jul 1 2022 11:17:40 for proteus by doxygen 1.8.20