1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 #include <ceed/types.h>
5 #include "utils.h"
6
7 /**
8 @brief Calculate CFL for 3D spatial dimensions
9
10 Defined as
11 @f[
12 \mathrm{CFL} = \Delta t \sqrt{\overline{g}_{jk} u_j u_k}
13 @f]
14
15 where \f$ \Delta t \f$ is the `timestep`, \f$ \overline{g}_{jk} \f$ is the element metric tensor `gij`, and \f$ u_j \f$ is the `velocity`.
16 For best results, `gij` should be scaled to account for the reference element domain (often [-1,1] for Gaussian quadrature).
17
18 @param[in] velocity Velocity
19 @param[in] timestep Timestep
20 @param[in] gij Element metric tensor (should be scaled)
21 @return The CFL value for the element
22 **/
CalculateCFL_3D(const CeedScalar velocity[3],CeedScalar timestep,const CeedScalar gij[3][3])23 CEED_QFUNCTION_HELPER CeedScalar CalculateCFL_3D(const CeedScalar velocity[3], CeedScalar timestep, const CeedScalar gij[3][3]) {
24 CeedScalar gij_uj[3] = {0.};
25
26 MatVec3(gij, velocity, CEED_NOTRANSPOSE, gij_uj);
27 return sqrt(Dot3(velocity, gij_uj)) * timestep;
28 }
29
30 /**
31 @brief Calculate CFL for 2D spatial dimensions
32
33 Defined as
34 @f[
35 \mathrm{CFL} = \Delta t \sqrt{\overline{g}_{jk} u_j u_k}
36 @f]
37
38 where \f$ \Delta t \f$ is the `timestep`, \f$ \overline{g}_{jk} \f$ is the element metric tensor `gij`, and \f$ u_j \f$ is the `velocity`.
39 For best results, `gij` should be scaled to account for the reference element domain (often [-1,1] for Gaussian quadrature).
40
41 @param[in] velocity Advection velocity
42 @param[in] timestep Timestep
43 @param[in] gij Element metric tensor (should be scaled)
44 @return The CFL value for the element
45 **/
CalculateCFL_2D(const CeedScalar velocity[2],CeedScalar timestep,const CeedScalar gij[2][2])46 CEED_QFUNCTION_HELPER CeedScalar CalculateCFL_2D(const CeedScalar velocity[2], CeedScalar timestep, const CeedScalar gij[2][2]) {
47 CeedScalar gij_uj[2] = {0.};
48
49 MatVec2(gij, velocity, CEED_NOTRANSPOSE, gij_uj);
50 return sqrt(Dot2(velocity, gij_uj)) * timestep;
51 }
52
53 /**
54 @brief Calculate Peclet number for 3D spatial dimensions
55
56 Defined as
57 @f[
58 \mathrm{Pe} = \frac{\sqrt{\overline{g}_{jk} u_j u_k}}{\kappa}
59 @f]
60
61 where \f$ \kappa \f$ is the `diffusion_coeff`, \f$ \overline{g}_{jk} \f$ is the element metric tensor `gij`, and \f$ u_j \f$ is the `velocity`.
62 For best results, `gij` should be scaled to account for the reference element domain (often [-1,1] for Gaussian quadrature).
63
64 @param[in] velocity Advection velocity
65 @param[in] diffusion_coeff Diffusion coefficient
66 @param[in] gij Element metric tensor (should be scaled)
67 @return The Peclet value for the element
68 **/
CalculatePe_3D(const CeedScalar velocity[3],CeedScalar diffusion_coeff,const CeedScalar gij[3][3])69 CEED_QFUNCTION_HELPER CeedScalar CalculatePe_3D(const CeedScalar velocity[3], CeedScalar diffusion_coeff, const CeedScalar gij[3][3]) {
70 CeedScalar gij_uj[3] = {0.}, gij_inv_mat[3][3] = {{0.}};
71
72 MatInv3(gij, gij_inv_mat, NULL);
73 MatVec3(gij_inv_mat, velocity, CEED_NOTRANSPOSE, gij_uj);
74 return sqrt(Dot3(velocity, gij_uj)) / diffusion_coeff;
75 }
76
77 /**
78 @brief Calculate Peclet number for 2D spatial dimensions
79
80 Defined as
81 @f[
82 \mathrm{Pe} = \frac{\sqrt{\overline{g}_{jk} u_j u_k}}{\kappa}
83 @f]
84
85 where \f$ \kappa \f$ is the `diffusion_coeff`, \f$ \overline{g}_{jk} \f$ is the element metric tensor `gij`, and \f$ u_j \f$ is the `velocity`.
86 For best results, `gij` should be scaled to account for the reference element domain (often [-1,1] for Gaussian quadrature).
87
88 @param[in] velocity Advection velocity
89 @param[in] diffusion_coeff Diffusion coefficient
90 @param[in] gij Element metric tensor (should be scaled)
91 @return The Peclet value for the element
92 **/
CalculatePe_2D(const CeedScalar velocity[2],CeedScalar diffusion_coeff,const CeedScalar gij[2][2])93 CEED_QFUNCTION_HELPER CeedScalar CalculatePe_2D(const CeedScalar velocity[2], CeedScalar diffusion_coeff, const CeedScalar gij[2][2]) {
94 CeedScalar gij_uj[2] = {0.}, gij_inv_mat[2][2] = {{0.}};
95
96 MatInv2(gij, gij_inv_mat, NULL);
97 MatVec2(gij_inv_mat, velocity, CEED_NOTRANSPOSE, gij_uj);
98 return sqrt(Dot2(velocity, gij_uj)) / diffusion_coeff;
99 }
100