// // Created by jundafeng on 3/23/22. // #ifndef DGTD_KERNELS_CU #define DGTD_KERNELS_CU #include "kernels.cuh" /* Basic kernel: one tetra per thread * Opt1: one block per tetra, one entry per thread, registers can be shuffled. * New kernels can be parallelized in a smarter way: * 1. volume integral can be parallelized in (12,12) * 2. surface integral can be parallelized in (6,6) * 3. coupling integral can be parallelized in (6,6) * 4. excitation can be parallelized in (12,1) */ //-----------------------CONSTANTS-------------------------------------------------------------------------------------------- __device__ constexpr int GAUSS_POINT_NUM[4] = {6, 9, 9, 9}; __device__ constexpr int TetPolyOrderDim[4] = {6, 12, 30, 45}; __device__ constexpr int FacePolyOrderDim[4] = {3, 6, 12, 18}; #define g6_a 0.816847572980459 #define g6_b (1.0 - g6_a) / 2.0 #define g6_c 0.108103018168070 #define g6_d (1.0 - g6_c) / 2.0 #define g6_W1 0.109951743655322 #define g6_W2 0.223381589678011 #define g9_a 0.437525248383384 #define g9_b 1.0 - 2.0 * g9_a #define g9_c 0.797112651860071 #define g9_d 0.165409927389841 #define g9_e 1.0 - g9_c - g9_d #define g9_W1 0.205950504760887 #define g9_W2 0.063691414286223 #if defined(DGTD_USE_CUDA) #include #include #include #include "cuda_utils.h" __device__ fp_t_ts g2d_6[6][4] = { g6_a, g6_b, g6_b, g6_W1, g6_b, g6_a, g6_b, g6_W1, g6_b, g6_b, g6_a, g6_W1, g6_c, g6_d, g6_d, g6_W2, g6_d, g6_c, g6_d, g6_W2, g6_d, g6_d, g6_c, g6_W2}; // added by Seung-Cheol Lee & Marinos N. Vouvakis 08/12/02 // Order = 5 __device__ fp_t_ts g2d_9[9][4] = { g9_b, g9_a, g9_a, g9_W1, g9_a, g9_b, g9_a, g9_W1, g9_a, g9_a, g9_b, g9_W1, g9_c, g9_d, g9_e, g9_W2, g9_c, g9_e, g9_d, g9_W2, g9_d, g9_c, g9_e, g9_W2, g9_d, g9_e, g9_c, g9_W2, g9_e, g9_c, g9_d, g9_W2, g9_e, g9_d, g9_c, g9_W2}; // __device__ const int FACE_BASIS_IJK[6][2] = { // {1, 2}, // {0, 2}, // {0, 1}, // {1, 2}, // {0, 2}, // {0, 1}, // }; __device__ const int FACE_BASIS_IJK[18][3] = { {1, 2, -1}, {0, 2, -1}, {0, 1, -1}, {1, 2, -1}, {0, 2, -1}, {0, 1, -1}, {0, 1, 2}, {0, 1, 2}, {1, 2, -1}, {0, 2, -1}, {0, 1, -1}, {0, 1, 2}, {0, 1, 2}, {0, 1, 2}, {0, 1, 2}, {0, 1, 2}, {0, 1, 2}, {0, 1, 2} }; __device__ const int FAC2TET[4][18] = { {5, 4, 3, 11, 10, 9, 12, 13, 25, 24, 23, 26, 30, 31, 32, 42, 43, 44}, {5, 2, 1, 11, 8, 7, 14, 15, 25, 22, 21, 27, 33, 34, 35, 42, 43, 44}, {4, 2, 0, 10, 8, 6, 16, 17, 24, 22, 20, 28, 36, 37, 38, 42, 43, 44}, {3, 1, 0, 9, 7, 6, 18, 19, 23, 21, 20, 29, 39, 40, 41, 42, 43, 44} }; #endif //---------------------------------------------------------------------------------------------------------------------------- #if defined(DGTD_USE_CUDA) || defined(DGTD_USE_CUDA_OPENCL) __global__ void makeNeighField(int* MappingArray, fp_t_ts* Field, fp_t_ts* auxiliar, int elements) { int idx = blockIdx.x * blockDim.x + threadIdx.x; if(idx < elements){ auxiliar[idx] = Field[MappingArray[idx]]; // printf("%d\n", MappingArray[idx]); } } __global__ void makeNeighFieldRegular(int* MappingArray, fp_t_ts* Field, fp_t_ts* auxiliar, int elements, int polyOrder){ int idx = blockIdx.x * blockDim.x + threadIdx.x; if(idx < NumOfFaces * FacePolyOrderDim[polyOrder] * elements){ int blockId = idx / FacePolyOrderDim[polyOrder]; int basisId = idx % FacePolyOrderDim[polyOrder]; int tetId = blockId / NumOfFaces; int faceId = blockId % NumOfFaces; auxiliar[(faceId * elements + tetId) * FacePolyOrderDim[polyOrder] + basisId] = Field[MappingArray[idx]]; // for(int i = 0; i < NumOfFaces; i++){ // auxiliar[idx + i * elements] = Field[MappingArray[idx + i * FacePolyOrderDim[polyOrder]]]; // } } } __global__ void makeNeighFieldRegular_custom(int* MappingArray, fp_t_ts* Field, fp_t_ts* auxiliar, int elements, int polyOrder, int numFaces){ int idx = blockIdx.x * blockDim.x + threadIdx.x; if(idx < numFaces * FacePolyOrderDim[polyOrder] * elements){ int blockId = idx / FacePolyOrderDim[polyOrder]; int basisId = idx % FacePolyOrderDim[polyOrder]; int tetId = blockId / numFaces; int faceId = blockId % numFaces; auxiliar[(faceId * elements + tetId) * FacePolyOrderDim[polyOrder] + basisId] = Field[MappingArray[idx]]; } } __global__ void addCouplingResults(fp_t_ts* Field, fp_t_ts* auxiliar, int* neighs, int* neighsOffset, int elements){ int tet = blockIdx.x * blockDim.y + threadIdx.y; if(tet < elements){ int basis = blockDim.x; fp_t_ts fieldAux = 0.0; int offset = (neighsOffset[tet] - neighsOffset[0]) * basis; for(int i = 0; i < neighs[tet]; i++){ fieldAux += auxiliar[offset + basis * i + threadIdx.x]; } Field[tet * basis + threadIdx.x] += fieldAux; } } __global__ void addCouplingResultsRegular(fp_t_ts* Field, fp_t_ts* auxiliar, int elements){ int tet = blockIdx.x * blockDim.y + threadIdx.y; if(tet < elements){ int basis = blockDim.x; fp_t_ts fieldAux = 0.0; for(int i = 0; i < NumOfFaces; i++){ fieldAux += auxiliar[tet * basis + elements * basis * i + threadIdx.x]; } Field[tet * basis + threadIdx.x] += fieldAux; } } __global__ void addCouplingResultsRegular_custom(fp_t_ts* Field, fp_t_ts* auxiliar, int elements, int numFaces){ int tet = blockIdx.x * blockDim.y + threadIdx.y; if(tet < elements) { int basis = blockDim.x; fp_t_ts fieldAux = 0.0; for(int i = 0; i < numFaces; i++) { fieldAux += auxiliar[tet * basis + elements * basis * i + threadIdx.x]; } Field[tet * basis + threadIdx.x] += fieldAux; } } __global__ void fill_ones_cols(fp_t_ts* __restrict__ M, int rows, int cols, int ld) { int r = blockIdx.x * blockDim.x + threadIdx.x; int c = blockIdx.y * blockDim.y + threadIdx.y; if (r < rows && c < cols) M[(size_t)c * ld + r] = 1.0f; } // 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 // // 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 // // 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 // __device__ fp_t_ts dotP3131(const fp_t_ts *v1, const fp_t_ts *v2){ fp_t_ts result; result = v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]; return result; } __device__ void crossP3(const fp_t_ts *v1, const fp_t_ts *v2, fp_t_ts *result){ result[0] = v1[1] * v2[2] - v1[2] * v2[1]; result[1] = -(v1[0] * v2[2] - v1[2] * v2[0]); result[2] = v1[0] * v2[1] - v1[1] * v2[0]; } __device__ void scale3(fp_t_ts *v1, const fp_t_ts factor){ v1[0] = v1[0] * factor; v1[1] = v1[1] * factor; v1[2] = v1[2] * factor; } __device__ void scale3(const fp_t_ts *v1, const fp_t_ts factor, fp_t_ts *v2){ v2[0] = v1[0] * factor; v2[1] = v1[1] * factor; v2[2] = v1[2] * factor; } __device__ void unit3(fp_t_ts *v1){ fp_t_ts work[3]; fp_t_ts mag = sqrt(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]); work[0] = v1[0] / mag; work[1] = v1[1] / mag; work[2] = v1[2] / mag; v1[0] = work[0]; v1[1] = work[1]; v1[2] = work[2]; } __device__ void plus3(fp_t_ts *v1, fp_t_ts *v2){ v1[0] = v1[0] + v2[0]; v1[1] = v1[1] + v2[1]; v1[2] = v1[2] + v2[2]; } __device__ void set3(const fp_t_ts *v1, fp_t_ts *v2){ v2[0] = v1[0]; v2[1] = v1[1]; v2[2] = v1[2]; } __device__ void plus3(const fp_t_ts *v1, const fp_t_ts *v2, fp_t_ts *result){ result[0] = v1[0] + v2[0]; result[1] = v1[1] + v2[1]; result[2] = v1[2] + v2[2]; } __device__ void plus33(const fp_t_ts *v1, const fp_t_ts *v2, const fp_t_ts *v3, fp_t_ts *result){ result[0] = v1[0] + v2[0] + v3[0]; result[1] = v1[1] + v2[1] + v3[1]; result[2] = v1[2] + v2[2] + v3[2]; } __device__ void minus3(const fp_t_ts *v1, const fp_t_ts *v2, fp_t_ts *result){ result[0] = v1[0] - v2[0]; result[1] = v1[1] - v2[1]; result[2] = v1[2] - v2[2]; } __device__ void reset3(fp_t_ts *result){ result[0] = 0.0; result[1] = 0.0; result[2] = 0.0; } __device__ void geometry3(const fp_t_ts *nd_coord, fp_t_ts lvtr[][3], fp_t_ts avtr[][3], fp_t_ts *vol){ int i, i1, i2, i3; fp_t_ts vv[3], tt[3]; for (i = 0; i < 3; i++){ minus3(&nd_coord[(i + 1) * 3], &nd_coord[0], lvtr[i]); } for (i = 0; i < 3; i++){ i1 = (i + 1) % 4; i2 = (i + 2) % 4; i3 = (i + 3) % 4; minus3(&nd_coord[i2 * 3], &nd_coord[i1 * 3], tt); minus3(&nd_coord[i3 * 3], &nd_coord[i1 * 3], vv); crossP3(tt, vv, avtr[i]); scale3(avtr[i], 0.5); minus3(&nd_coord[i * 3], &nd_coord[i1 * 3], tt); if (dotP3131(avtr[i], tt) < 0.0) scale3(avtr[i], -1); ; } reset3(avtr[3]); plus33(avtr[0], avtr[1], avtr[2], avtr[3]); scale3(avtr[3], -1); (*vol) = dotP3131(avtr[2], lvtr[1]) / 3.0; } __device__ fp_t_ts magnitude3(const fp_t_ts *v1){ return sqrt(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]); } __device__ void tet_center(const fp_t_ts *tet_nd_coord, fp_t_ts *tet_center){ tet_center[0] = (tet_nd_coord[0 * 3 + 0] + tet_nd_coord[1 * 3 + 0] + tet_nd_coord[2 * 3 + 0] + tet_nd_coord[3 * 3 + 0]) * 0.25; tet_center[1] = (tet_nd_coord[0 * 3 + 1] + tet_nd_coord[1 * 3 + 1] + tet_nd_coord[2 * 3 + 1] + tet_nd_coord[3 * 3 + 1]) * 0.25; tet_center[2] = (tet_nd_coord[0 * 3 + 2] + tet_nd_coord[1 * 3 + 2] + tet_nd_coord[2 * 3 + 2] + tet_nd_coord[3 * 3 + 2]) * 0.25; } // have to use 2d arrays in the device functions, using 1d array with offsets will result in erroneous values __device__ void geometry2(const fp_t_ts *nd_coord, const fp_t_ts *tet0_nd_coord, fp_t_ts *normal, fp_t_ts nvtr[][3], fp_t_ts *area){ fp_t_ts t0[3], t1[3], tv[3], tet0center[3]; minus3(&nd_coord[2 * 3], &nd_coord[1 * 3], t0); minus3(&nd_coord[0 * 3], &nd_coord[2 * 3], t1); crossP3(t0, t1, normal); (*area) = 0.5 * magnitude3(normal); scale3(normal, 1 / magnitude3(normal)); crossP3(normal, t0, nvtr[0]); crossP3(normal, t1, nvtr[1]); if (dotP3131(nvtr[0], t1) < 0.0) scale3(nvtr[0], -1); if (dotP3131(nvtr[1], t0) > 0.0) scale3(nvtr[1], -1); tet_center(tet0_nd_coord, tet0center); minus3(tet0center, &nd_coord[0 * 3], tv); if (dotP3131(tv, normal) > 0.0) scale3(normal, -1); } __device__ void geometry2(const fp_t_ts *nd_coord, fp_t_ts nvtr[][3], fp_t_ts *area){ fp_t_ts t0[3], t1[3], normal[3]; minus3(&nd_coord[2 * 3], &nd_coord[1 * 3], t0); minus3(&nd_coord[0 * 3], &nd_coord[2 * 3], t1); crossP3(t0, t1, normal); (*area) = 0.5 * magnitude3(normal); scale3(normal, 1 / magnitude3(normal)); crossP3(normal, t0, nvtr[0]); crossP3(normal, t1, nvtr[1]); if (dotP3131(nvtr[0], t1) < 0.0) scale3(nvtr[0], -1); if (dotP3131(nvtr[1], t0) > 0.0) scale3(nvtr[1], -1); } __device__ void dotP3331(const fp_t_ts *tensor, const fp_t_ts *v1, fp_t_ts *result){ result[0] = v1[0] * tensor[0] + v1[1] * tensor[1] + v1[2] * tensor[2]; result[1] = v1[0] * tensor[3] + v1[1] * tensor[4] + v1[2] * tensor[5]; result[2] = v1[0] * tensor[6] + v1[1] * tensor[7] + v1[2] * tensor[8]; } __device__ void makeCoeff33(const fp_t_ts avtr[][3], const fp_t_ts *epsr, fp_t_ts coeffA[][3]){ int i, j; fp_t_ts temp[3]; for (i = 0; i < 3; i++){ for (j = 0; j < 3; j++){ dotP3331(epsr, avtr[j], temp); coeffA[i][j] = dotP3131(temp, avtr[i]); } } } __device__ void quadObtention(int np, int id, fp_t_ts &zeta0, fp_t_ts &zeta1, fp_t_ts &zeta2, fp_t_ts &weight){ if (id < 0 || id >= np){ printf("Error in quadObtention"); return; } switch (np){ case 6: zeta0 = g2d_6[id][0]; zeta1 = g2d_6[id][1]; zeta2 = g2d_6[id][2]; weight = g2d_6[id][3]; break; case 9: zeta0 = g2d_9[id][0]; zeta1 = g2d_9[id][1]; zeta2 = g2d_9[id][2]; weight = g2d_9[id][3]; break; default: break; } } __device__ void faceBasisW(const fp_t_ts *node_coords, const fp_t_ts *zeta, const int p, fp_t_ts *W){ fp_t_ts gradZeta[3][3]; fp_t_ts area; geometry2(node_coords, gradZeta, &area); scale3(gradZeta[0], 0.5 / area); scale3(gradZeta[1], 0.5 / area); gradZeta[2][0] = 0; gradZeta[2][1] = 0; gradZeta[2][2] = 0; minus3(gradZeta[2], gradZeta[0], gradZeta[2]); minus3(gradZeta[2], gradZeta[1], gradZeta[2]); int i = FACE_BASIS_IJK[p][0]; int j = FACE_BASIS_IJK[p][1]; int k = FACE_BASIS_IJK[p][2]; switch (p){ case 0: scale3(gradZeta[j], zeta[i]); scale3(gradZeta[i], zeta[j]); minus3(gradZeta[j], gradZeta[i], W); break; case 1: scale3(gradZeta[j], zeta[i]); scale3(gradZeta[i], zeta[j]); minus3(gradZeta[j], gradZeta[i], W); break; case 2: scale3(gradZeta[j], zeta[i]); scale3(gradZeta[i], zeta[j]); minus3(gradZeta[j], gradZeta[i], W); break; case 3: scale3(gradZeta[j], zeta[i]); scale3(gradZeta[i], zeta[j]); plus3(gradZeta[j], gradZeta[i], W); scale3(W, 4.); break; case 4: scale3(gradZeta[j], zeta[i]); scale3(gradZeta[i], zeta[j]); plus3(gradZeta[j], gradZeta[i], W); scale3(W, 4.); break; case 5: scale3(gradZeta[j], zeta[i]); scale3(gradZeta[i], zeta[j]); plus3(gradZeta[j], gradZeta[i], W); scale3(W, 4.); break; case 6: scale3(gradZeta[k], zeta[i]); scale3(gradZeta[i], zeta[k]); minus3(gradZeta[k], gradZeta[i], W); scale3(W, 4 * zeta[j]); break; case 7: scale3(gradZeta[j], zeta[i]); scale3(gradZeta[i], zeta[j]); minus3(gradZeta[j], gradZeta[i], W); scale3(W, 4. * zeta[k]); break; case 8: scale3(gradZeta[j], zeta[i] * zeta[i] - 2. * zeta[i] * zeta[j]); scale3(gradZeta[i], zeta[j] * zeta[j] - 2. * zeta[i] * zeta[j]); minus3(gradZeta[j], gradZeta[i], W); scale3(W, 4.); break; case 9: scale3(gradZeta[j], zeta[i] * zeta[i] - 2. * zeta[i] * zeta[j]); scale3(gradZeta[i], zeta[j] * zeta[j] - 2. * zeta[i] * zeta[j]); minus3(gradZeta[j], gradZeta[i], W); scale3(W, 4.); break; case 10: scale3(gradZeta[j], zeta[i] * zeta[i] - 2. * zeta[i] * zeta[j]); scale3(gradZeta[i], zeta[j] * zeta[j] - 2. * zeta[i] * zeta[j]); minus3(gradZeta[j], gradZeta[i], W); scale3(W, 4.); break; case 11: scale3(gradZeta[k], zeta[i] * zeta[j]); scale3(gradZeta[j], zeta[i] * zeta[k]); scale3(gradZeta[i], zeta[k] * zeta[j]); plus33(gradZeta[k], gradZeta[j], gradZeta[i], W); scale3(W, 8.); break; case 12: scale3(gradZeta[k], zeta[j]); scale3(gradZeta[j], zeta[k]); minus3(gradZeta[k], gradZeta[j], W); scale3(W, 4. * zeta[i] * zeta[i]); break; case 13: scale3(gradZeta[k], zeta[i]); scale3(gradZeta[i], zeta[k]); minus3(gradZeta[k], gradZeta[i], W); scale3(W, 4. * zeta[j] * zeta[j]); break; case 14: scale3(gradZeta[j], zeta[i]); scale3(gradZeta[i], zeta[j]); minus3(gradZeta[j], gradZeta[i], W); scale3(W, 4. * zeta[k] * zeta[k]); break; case 15: W[0] = W[1] = W[2] = 0; break; case 16: W[0] = W[1] = W[2] = 0; break; case 17: W[0] = W[1] = W[2] = 0; break; default: W[0] = W[1] = W[2] = 0; printf("Check faceBasisW: you are asking for somthing out of scope"); break; } } __device__ void getCartesianCoord(const fp_t_ts *node_coords, const fp_t_ts *zeta, fp_t_ts *pnt){ pnt[0] = node_coords[0 * 3 + 0] * zeta[0] + node_coords[1 * 3 + 0] * zeta[1] + node_coords[2 * 3 + 0] * zeta[2]; pnt[1] = node_coords[0 * 3 + 1] * zeta[0] + node_coords[1 * 3 + 1] * zeta[1] + node_coords[2 * 3 + 1] * zeta[2]; pnt[2] = node_coords[0 * 3 + 2] * zeta[0] + node_coords[1 * 3 + 2] * zeta[1] + node_coords[2 * 3 + 2] * zeta[2]; } __device__ void PlaneWaveFields(int ExcitationFlag, fp_t_ts t, fp_t_ts to, fp_t_ts tau, const fp_t_ts *r, const fp_t_ts *ro, fp_t_ts freq_m, fp_t_ts *kvtr, fp_t_ts eta, fp_t_ts Emagnitude, fp_t_ts &IncidExcit_E, fp_t_ts &IncidExcit_H) { fp_t_ts omega = 2.0 * Pi * freq_m * MEGA; fp_t_ts Exponent; fp_t_ts SinModul; fp_t_ts Neuman; fp_t_ts work[3]; // printf("%d", ExcitationFlag); switch (ExcitationFlag){ case 0: minus3(r, ro, work); Exponent = t - to - dotP3131(kvtr, work) / Vo; // Time Harmonic if (Exponent >= 0.0){ SinModul = cos(omega * Exponent); IncidExcit_E = Emagnitude * SinModul; IncidExcit_H = (Emagnitude / eta) * SinModul; }else{ IncidExcit_E = 0.0; IncidExcit_H = 0.0; } break; case 1: // Gauss Pulse minus3(r, ro, work); Exponent = t - to - dotP3131(kvtr, work) / Vo; SinModul = cos(omega * Exponent); IncidExcit_E = Emagnitude * SinModul * exp(-(Exponent * Exponent) / (tau * tau)); IncidExcit_H = (Emagnitude / eta) * SinModul * exp(-(Exponent * Exponent) / (tau * tau)); break; case 2: // Neuman Pulse minus3(r, ro, work); Exponent = t - to - dotP3131(kvtr, work) / Vo; Neuman = (2.0 * Exponent) / (tau * tau); IncidExcit_E = (Emagnitude * Neuman) * exp(-(Exponent * Exponent) / (tau * tau)); IncidExcit_H = ((Emagnitude / eta) * Neuman) * exp(-(Exponent * Exponent) / (tau * tau)); break; case 3: { fp_t_ts work[3]; minus3(r, ro, work); // work = r - ro Exponent = t - dotP3131(kvtr, work) / Vo; if (Exponent >= 0.0 && Exponent <= to) { // Chirp parameters fp_t_ts f_end = freq_m * MEGA; // Center frequency in Hz fp_t_ts B = tau * MEGA; // Bandwidth in Hz fp_t_ts f0 = f_end - B; // Start frequency fp_t_ts f1 = f_end; // End frequency fp_t_ts Tchirp = to; // Chirp duration // Chirp phase fp_t_ts chirpArg = 2.0 * Pi * f0 * Exponent + Pi * (f1 - f0) / Tchirp * Exponent * Exponent; // Final tapered chirp signal fp_t_ts chirpSignal = sin(chirpArg); IncidExcit_E = Emagnitude * chirpSignal; IncidExcit_H = (Emagnitude / eta) * chirpSignal; } else { IncidExcit_E = 0.0; IncidExcit_H = 0.0; } break; } default: break; } } __global__ void CalculatePadeFreq(fp_t* a_k, fp_t* b_k, int M_global, int N, int* constantM, cuDoubleComplex* Hf){ int m = blockIdx.x * blockDim.x + threadIdx.x; if(m < M_global){ fp_t sum_real_a = 0.0; fp_t sum_imag_a = 0.0; fp_t sum_real_b = 0.0; fp_t sum_imag_b = 0.0; for (int n = 0; n < N; n++) { fp_t angle = -2.0f * Pi * constantM[m] * n / M_global; fp_t cos_val, sin_val; sincos(angle, &sin_val, &cos_val); sum_real_a += a_k[n] * cos_val; sum_imag_a += a_k[n] * sin_val; sum_real_b += b_k[n] * cos_val; sum_imag_b += b_k[n] * sin_val; } Hf[m] = cuCdiv(make_cuDoubleComplex(sum_real_a, sum_imag_a), make_cuDoubleComplex(sum_real_b, sum_imag_b)); } return; } // === Utility Functions for Double Precision === __device__ __forceinline__ fp_t_ts cleanSmall(fp_t_ts val, fp_t_ts eps = 1e-12) { return fabs(val) < eps ? 0.0 : val; } __device__ __forceinline__ void cleanVec3(fp_t_ts* v, fp_t_ts eps = 1e-12) { for (int i = 0; i < 3; ++i) v[i] = cleanSmall(v[i], eps); } // -------------------------------------------------------------------------------------- // Plane Wave Exciation Update // -------------------------------------------------------------------------------------- __device__ void E_updateinc_PW(const fp_t_ts t, const fp_t_ts to, const fp_t_ts tau, const fp_t_ts freq_m, const fp_t_ts eta, const fp_t_ts Emagnitude, const fp_t_ts *ro, const fp_t_ts theta, const fp_t_ts phi, const fp_t_ts *Epol, const int ExcitationFlag, const fp_t_ts *r, const fp_t_ts *normalvtr, fp_t_ts *nxHinc, fp_t_ts *nxnxEinc) { fp_t_ts IncidExcit_E, IncidExcit_H; fp_t_ts theta_in_rad = theta * (fp_t_ts)Pi / (fp_t_ts)180.0; fp_t_ts phi_in_rad = phi * (fp_t_ts)Pi / (fp_t_ts)180.0; fp_t_ts kvtr[3], Hpol[3], work[3]; kvtr[0] = sin(theta_in_rad) * cos(phi_in_rad); kvtr[1] = sin(theta_in_rad) * sin(phi_in_rad); kvtr[2] = cos(theta_in_rad); // Fields in Time Domain PlaneWaveFields(ExcitationFlag, t, to, tau, r, ro, freq_m, kvtr, eta, Emagnitude, IncidExcit_E, IncidExcit_H); // EField crossP3(normalvtr, Epol, work); crossP3(normalvtr, work, nxnxEinc); scale3(nxnxEinc, IncidExcit_E); // HField crossP3(kvtr, Epol, work); crossP3(normalvtr, work, Hpol); scale3(Hpol, IncidExcit_H, nxHinc); /* if (r[2] == 0.0 && r[1] < 1e-4 && r[1] > 0.0 && r[0] < 5e-4 && r[0] > 0.0) { printf( "r=(%e,%e,%e) | nxnxE =(%e,%e,%e) | nxH =(%e,%e,%e) \n", (double)r[0], (double)r[1], (double)r[2], (double)nxnxEinc[0], (double)nxnxEinc[1], (double)nxnxEinc[2], (double)nxHinc[0], (double)nxHinc[1], (double)nxHinc[2] ); } */ for (int i = 0; i < 3; ++i) { if (fabs(nxHinc[i]) < (fp_t_ts)1e-12) nxHinc[i] = 0.0; if (fabs(nxnxEinc[i]) < (fp_t_ts)1e-12) nxnxEinc[i] = 0.0; } } __device__ void H_updateinc_PW(const fp_t_ts t, const fp_t_ts to, const fp_t_ts tau, const fp_t_ts freq_m, const fp_t_ts eta, const fp_t_ts Emagnitude, const fp_t_ts *ro, const fp_t_ts theta, const fp_t_ts phi, const fp_t_ts *Epol, const int ExcitationFlag, const fp_t_ts *r, const fp_t_ts *normalvtr, fp_t_ts *nxEinc, fp_t_ts *nxnxHinc) { fp_t_ts IncidExcit_E, IncidExcit_H; fp_t_ts theta_in_rad = theta * (fp_t_ts)Pi / (fp_t_ts)180.0; fp_t_ts phi_in_rad = phi * (fp_t_ts)Pi / (fp_t_ts)180.0; fp_t_ts kvtr[3], Hpol[3], work[3]; kvtr[0] = sin(theta_in_rad) * cos(phi_in_rad); kvtr[1] = sin(theta_in_rad) * sin(phi_in_rad); kvtr[2] = cos(theta_in_rad); PlaneWaveFields(ExcitationFlag, t, to, tau, r, ro, freq_m, kvtr, eta, Emagnitude, IncidExcit_E, IncidExcit_H); // EField crossP3(normalvtr, Epol, work); scale3(work, IncidExcit_E, nxEinc); // HField crossP3(kvtr, Epol, Hpol); crossP3(normalvtr, Hpol, work); crossP3(normalvtr, work, nxnxHinc); scale3(nxnxHinc, IncidExcit_H); for (int i = 0; i < 3; ++i) { if (fabs(nxEinc[i]) < (fp_t_ts)1e-12) nxEinc[i] = 0.0; if (fabs(nxnxHinc[i]) < (fp_t_ts)1e-12) nxnxHinc[i] = 0.0; } } __global__ void addExcitationH_PML( int* ExcitationFacesCnt, int* ExcitationFacesOffset, int* ExcitationFacesNum, int excitationCNT, int excitationPMLCNT, int8_t* map, ExcitationProp exci_prop, int PolyFlag, fp_t_ts dt_muo, fp_t_ts t, fp_t_ts* nd_coords_tet, fp_t_ts* nd_coords_face, fp_t_ts* Z_face, fp_t_ts* InvMat, fp_t_ts* Field) { int tet_id = blockIdx.x * blockDim.x + threadIdx.x; if (tet_id < excitationCNT) { //-----TODO: Make this part dynamic---------------------------- const int TetPolyOrderDim2 = TetPolyOrderDim[3]; const int FacePolyOrderDim2 = FacePolyOrderDim[3]; fp_t_ts nxe_inc[TetPolyOrderDim2] = {0}; fp_t_ts nxnxh_inc[TetPolyOrderDim2] = {0}; //------------------------------------------------------------- fp_t_ts Z = 0.0; for (int i = 0; i < ExcitationFacesCnt[tet_id]; i++){ int face_id = ExcitationFacesOffset[tet_id] + i; Z = Z_face[face_id]; fp_t_ts nvtr[2][3]; fp_t_ts normvtr[3] = {0}; fp_t_ts area = 0; fp_t_ts nxentry; fp_t_ts nxnxentry; fp_t_ts z[3], w[3], r[3]; fp_t_ts wgt; //-----------------TODO: Make this part dynamic-------------------------------------------- fp_t_ts nxelem[FacePolyOrderDim2] = {0}; fp_t_ts nxnxelem[FacePolyOrderDim2] = {0}; //----------------------------------------------------------------------------------------- fp_t_ts nxE_inc[3]; fp_t_ts nxnxH_inc[3]; geometry2(&nd_coords_face[3 * 3 * face_id], &nd_coords_tet[4 * 3 * tet_id], normvtr, nvtr, &area); for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++){ nxentry = 0; nxnxentry = 0; for (int k = 0; k < GAUSS_POINT_NUM[PolyFlag]; k++){ quadObtention(GAUSS_POINT_NUM[PolyFlag], k, z[0], z[1], z[2], wgt); faceBasisW(&nd_coords_face[3 * 3 * face_id], z, j, w); getCartesianCoord(&nd_coords_face[3 * 3 * face_id], z, r); H_updateinc_PW(t, exci_prop.to, exci_prop.tau, exci_prop.freq_m, Z, exci_prop.Emagnitude, exci_prop.ro, exci_prop.theta, exci_prop.phi, exci_prop.Epol, exci_prop.ExcitationFlag, r, normvtr, nxE_inc, nxnxH_inc); if (!(isnan(nxE_inc[0]))){ nxentry += area * wgt * dotP3131(w, nxE_inc); nxnxentry += area * wgt * dotP3131(w, nxnxH_inc); } } //nxelem[j] += nxentry; //nxnxelem[j] += nxnxentry; nxelem[j] += cleanSmall(nxentry); nxnxelem[j] += cleanSmall(nxnxentry); } for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++){ // printf("%d\n", ExcitationFacesNum[face_id]); nxe_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxelem[j]; nxnxh_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxnxelem[j]; } } for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++){ nxe_inc[i] *= -dt_muo; //nxe_inc[i] += (-Z * dt_muo) * nxnxh_inc[i]; nxe_inc[i] *= map[tet_id * TetPolyOrderDim[PolyFlag] + i]; } fp_t_ts sum; for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++){ sum = 0.0; for (int j = 0; j < TetPolyOrderDim[PolyFlag]; j++){ sum += InvMat[tet_id * (TetPolyOrderDim[PolyFlag] * TetPolyOrderDim[PolyFlag]) + j * TetPolyOrderDim[PolyFlag] + i] * nxe_inc[j] * map[tet_id * TetPolyOrderDim[PolyFlag] + i]; } sum = cleanSmall(sum); if (tet_id < excitationPMLCNT) { Field[tet_id * TetPolyOrderDim[PolyFlag] + i] -= 0.5 * sum; // Scattered Region } else { Field[tet_id * TetPolyOrderDim[PolyFlag] + i] += 0.5 * sum; // Total Region } } } } __global__ void addExcitationE_PML( int* ExcitationFacesCnt, int* ExcitationFacesOffset, int* ExcitationFacesNum, int excitationCNT, int excitationPMLCNT, int8_t* map, ExcitationProp exci_prop, int PolyFlag, fp_t_ts dt_epsilon, fp_t_ts t, fp_t_ts* nd_coords_tet, fp_t_ts* nd_coords_face, fp_t_ts* Z_face, fp_t_ts* InvMat, fp_t_ts* Field){ int tet_id = blockIdx.x * blockDim.x + threadIdx.x; if (tet_id < excitationCNT){ //-----TODO: Make this part dynamic---------------------------- const int TetPolyOrderDim2 = TetPolyOrderDim[3]; const int FacePolyOrderDim2 = FacePolyOrderDim[3]; fp_t_ts nxh_inc[TetPolyOrderDim2] = {0}; fp_t_ts nxnxe_inc[TetPolyOrderDim2] = {0}; //------------------------------------------------------------- fp_t_ts Y = 0.0; for (int i = 0; i < ExcitationFacesCnt[tet_id]; i++){ int face_id = ExcitationFacesOffset[tet_id] + i; Y = 1 / Z_face[face_id]; fp_t_ts nvtr[2][3]; fp_t_ts normvtr[3] = {0}; fp_t_ts area = 0; fp_t_ts nxentry; fp_t_ts nxnxentry; fp_t_ts z[3], w[3], r[3]; fp_t_ts wgt; //-----------------TODO: Make this part dynamic-------------------------------------------- fp_t_ts nxelem[FacePolyOrderDim2] = {0}; fp_t_ts nxnxelem[FacePolyOrderDim2] = {0}; //----------------------------------------------------------------------------------------- fp_t_ts nxH_inc[3]; fp_t_ts nxnxE_inc[3]; geometry2(&nd_coords_face[3 * 3 * face_id], &nd_coords_tet[4 * 3 * tet_id], normvtr, nvtr, &area); for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++){ nxentry = 0; nxnxentry = 0; for (int k = 0; k < GAUSS_POINT_NUM[PolyFlag]; k++){ quadObtention(GAUSS_POINT_NUM[PolyFlag], k, z[0], z[1], z[2], wgt); faceBasisW(&nd_coords_face[3 * 3 * face_id], z, j, w); getCartesianCoord(&nd_coords_face[3 * 3 * face_id], z, r); E_updateinc_PW(t, exci_prop.to, exci_prop.tau, exci_prop.freq_m, Z_face[face_id], exci_prop.Emagnitude, exci_prop.ro, exci_prop.theta, exci_prop.phi, exci_prop.Epol, exci_prop.ExcitationFlag, r, normvtr, nxH_inc, nxnxE_inc); if (!(isnan(nxH_inc[0]))){ nxentry += area * wgt * dotP3131(w, nxH_inc); nxnxentry += area * wgt * dotP3131(w, nxnxE_inc); } } //nxelem[j] += nxentry; //nxnxelem[j] += nxnxentry; nxelem[j] += cleanSmall(nxentry); nxnxelem[j] += cleanSmall(nxnxentry); } for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++){ // printf("%d\n", ExcitationFacesNum[face_id]); nxh_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxelem[j]; nxnxe_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxnxelem[j]; } } for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++){ nxh_inc[i] *= dt_epsilon; //nxh_inc[i] += (-Y * dt_epsilon) * nxnxe_inc[i]; nxh_inc[i] *= map[tet_id * TetPolyOrderDim[PolyFlag] + i]; } fp_t_ts sum; for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++){ sum = 0.0; for (int j = 0; j < TetPolyOrderDim[PolyFlag]; j++) { sum += InvMat[tet_id * (TetPolyOrderDim[PolyFlag] * TetPolyOrderDim[PolyFlag]) + j * TetPolyOrderDim[PolyFlag] + i] * nxh_inc[j] * map[tet_id * TetPolyOrderDim[PolyFlag] + i]; } sum = cleanSmall(sum); if (tet_id < excitationPMLCNT) { Field[tet_id * TetPolyOrderDim[PolyFlag] + i] -= 0.5 * sum; // Scattered Region } else { Field[tet_id * TetPolyOrderDim[PolyFlag] + i] += 0.5 * sum; // Total Region } } } } __global__ void addExcitationE( int* ExcitationFacesCnt, int* ExcitationFacesOffset, int* ExcitationFacesNum, int excitationCNT, int8_t* map, ExcitationProp exci_prop, int PolyFlag, fp_t_ts dt_epsilon, fp_t_ts t, fp_t_ts* nd_coords_tet, fp_t_ts* nd_coords_face, fp_t_ts* Z_face, fp_t_ts* InvMat, fp_t_ts* Field) { int tet_id = blockIdx.x * blockDim.x + threadIdx.x; if (tet_id < excitationCNT) { //-----TODO: Make this part dynamic---------------------------- const int TetPolyOrderDim2 = TetPolyOrderDim[3]; const int FacePolyOrderDim2 = FacePolyOrderDim[3]; fp_t_ts nxh_inc[TetPolyOrderDim2] = {0}; fp_t_ts nxnxe_inc[TetPolyOrderDim2] = {0}; //------------------------------------------------------------- fp_t_ts Y = 0.0; for (int i = 0; i < ExcitationFacesCnt[tet_id]; i++) { int face_id = ExcitationFacesOffset[tet_id] + i; Y = 1 / Z_face[face_id]; fp_t_ts nvtr[2][3]; fp_t_ts normvtr[3] = {0}; fp_t_ts area = 0; fp_t_ts nxentry; fp_t_ts nxnxentry; fp_t_ts z[3], w[3], r[3]; fp_t_ts wgt; //-----------------TODO: Make this part dynamic-------------------------------------------- fp_t_ts nxelem[FacePolyOrderDim2] = {0}; fp_t_ts nxnxelem[FacePolyOrderDim2] = {0}; //----------------------------------------------------------------------------------------- fp_t_ts nxH_inc[3]; fp_t_ts nxnxE_inc[3]; geometry2(&nd_coords_face[3 * 3 * face_id], &nd_coords_tet[4 * 3 * tet_id], normvtr, nvtr, &area); for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++) { nxentry = 0; nxnxentry = 0; for (int k = 0; k < GAUSS_POINT_NUM[PolyFlag]; k++) { quadObtention(GAUSS_POINT_NUM[PolyFlag], k, z[0], z[1], z[2], wgt); faceBasisW(&nd_coords_face[3 * 3 * face_id], z, j, w); getCartesianCoord(&nd_coords_face[3 * 3 * face_id], z, r); E_updateinc_PW(t, exci_prop.to, exci_prop.tau, exci_prop.freq_m, Z_face[face_id], exci_prop.Emagnitude, exci_prop.ro, exci_prop.theta, exci_prop.phi,exci_prop.Epol, exci_prop.ExcitationFlag,r,normvtr,nxH_inc,nxnxE_inc); if (!(isnan(nxH_inc[0]))) { nxentry += area * wgt * dotP3131(w, nxH_inc); nxnxentry += area * wgt * dotP3131(w, nxnxE_inc); } } //printf("face_id=%d j=%d nxentry=%.6e nxnxentry=%.6e\n", face_id, j, (double)nxentry, (double)nxnxentry); nxelem[j] += nxentry; nxnxelem[j] += nxnxentry; } for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++){ //printf("%d\n", ExcitationFacesNum[face_id]); nxh_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxelem[j]; nxnxe_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxnxelem[j]; //printf("tet_id=%d face_id=%d j=%d TetNode=%d nxh_inc=%.6e nxnxe_inc=%.6e\n", tet_id, face_id, j, FAC2TET[ExcitationFacesNum[face_id]][j], (double)nxh_inc[FAC2TET[ExcitationFacesNum[face_id]][j]], (double)nxnxe_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] ); } } for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++){ nxh_inc[i] *= dt_epsilon; //nxh_inc[i] += (-Y * dt_epsilon) * nxnxe_inc[i]; nxh_inc[i] += (-Y * dt_epsilon) * nxnxe_inc[i]; nxh_inc[i] *= map[tet_id * TetPolyOrderDim[PolyFlag] + i]; //printf("tet_id=%d i=%d nxh_inc=%.6e nxnxe_inc=%6e Y=%6e dt_epsilon=%6e \n", tet_id, i, (double)nxh_inc[i], (double)nxnxe_inc[i], (double)Y, (double)dt_epsilon ); } fp_t_ts sum; for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++){ sum = 0.0; for (int j = 0; j < TetPolyOrderDim[PolyFlag]; j++){ sum += InvMat[tet_id * (TetPolyOrderDim[PolyFlag] * TetPolyOrderDim[PolyFlag]) + j * TetPolyOrderDim[PolyFlag] + i] * nxh_inc[j] * map[tet_id * TetPolyOrderDim[PolyFlag] + i]; //printf("InvMat[%d][%d]=%.6e | nxh_inc[%d]=%.6e | map=%d\n", j, i, (double)InvMat[tet_id * (TetPolyOrderDim[PolyFlag] * TetPolyOrderDim[PolyFlag]) + j * TetPolyOrderDim[PolyFlag] + i], j, (double)nxh_inc[j], map[tet_id * TetPolyOrderDim[PolyFlag] + i]); } //printf("E-field update: tet_id=%d i=%d sum=%.6e\n", tet_id, i, (double)sum); Field[tet_id * TetPolyOrderDim[PolyFlag] + i] += sum; } } } __global__ void addExcitationH( int* ExcitationFacesCnt, int* ExcitationFacesOffset, int* ExcitationFacesNum, int excitationCNT, int8_t* map, ExcitationProp exci_prop, int PolyFlag, fp_t_ts dt_muo, fp_t_ts t, fp_t_ts* nd_coords_tet, fp_t_ts* nd_coords_face, fp_t_ts* Z_face, fp_t_ts* InvMat, fp_t_ts* Field) { int tet_id = blockIdx.x * blockDim.x + threadIdx.x; if (tet_id < excitationCNT) { //printf("ExcitationFacesOffset = %d\n",ExcitationFacesOffset[tet_id]); //-----TODO: Make this part dynamic---------------------------- const int TetPolyOrderDim2 = TetPolyOrderDim[3]; const int FacePolyOrderDim2 = FacePolyOrderDim[3]; fp_t_ts nxe_inc[TetPolyOrderDim2] = {0}; fp_t_ts nxnxh_inc[TetPolyOrderDim2] = {0}; //------------------------------------------------------------- fp_t_ts Z = 0.0; for (int i = 0; i < ExcitationFacesCnt[tet_id]; i++){ int face_id = ExcitationFacesOffset[tet_id] + i; Z = Z_face[face_id]; fp_t_ts nvtr[2][3]; fp_t_ts normvtr[3] = {0}; fp_t_ts area = 0; fp_t_ts nxentry; fp_t_ts nxnxentry; fp_t_ts z[3], w[3], r[3]; fp_t_ts wgt; //-----------------TODO: Make this part dynamic-------------------------------------------- fp_t_ts nxelem[FacePolyOrderDim2] = {0}; fp_t_ts nxnxelem[FacePolyOrderDim2] = {0}; //----------------------------------------------------------------------------------------- fp_t_ts nxE_inc[3]; fp_t_ts nxnxH_inc[3]; geometry2(&nd_coords_face[3 * 3 * face_id], &nd_coords_tet[4 * 3 * tet_id], normvtr, nvtr, &area); for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++) { nxentry = 0; nxnxentry = 0; for (int k = 0; k < GAUSS_POINT_NUM[PolyFlag]; k++) { quadObtention(GAUSS_POINT_NUM[PolyFlag], k, z[0], z[1], z[2], wgt); faceBasisW(&nd_coords_face[3 * 3 * face_id], z, j, w); getCartesianCoord(&nd_coords_face[3 * 3 * face_id], z, r); H_updateinc_PW(t, exci_prop.to, exci_prop.tau, exci_prop.freq_m, Z, exci_prop.Emagnitude, exci_prop.ro, exci_prop.theta, exci_prop.phi,exci_prop.Epol, exci_prop.ExcitationFlag,r,normvtr,nxE_inc,nxnxH_inc); if (!(isnan(nxE_inc[0]))) { nxentry += area * wgt * dotP3131(w, nxE_inc); nxnxentry += area * wgt * dotP3131(w, nxnxH_inc); /* if (tet_id == 6) { printf("area = %0.6e | wgt = %0.6e | w = %0.6e,%0.6e,%0.6e | nxE_inc = %0.6e,%0.6e,%0.6e | nxnxH_inc = %0.6e,%0.6e,%0.6e \n", area,wgt, w[0], w[1], w[2],nxE_inc[0],nxE_inc[1],nxE_inc[2],nxnxH_inc[0],nxnxH_inc[1],nxnxH_inc[2]); } */ } } nxelem[j] += nxentry; nxnxelem[j] += nxnxentry; } for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++) { // printf("%d\n", ExcitationFacesNum[face_id]); nxe_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxelem[j]; nxnxh_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxnxelem[j]; } } for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++) { /* if (tet_id == 3) { printf("tet=%d i=%d | nxe_inc=%.6e | dt_muo=%.6e Z=%.6e | nxnxh=%.6e | map[%d]=%d\n", tet_id, i, (double)nxe_inc[i], (double)dt_muo, (double)Z, (double)nxnxh_inc[i], map[tet_id * TetPolyOrderDim[PolyFlag] + i]); } */ nxe_inc[i] *= -dt_muo; nxe_inc[i] += (-Z * dt_muo) * nxnxh_inc[i]; nxe_inc[i] *= map[tet_id * TetPolyOrderDim[PolyFlag] + i]; } fp_t_ts sum; for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++) { sum = 0.0; for (int j = 0; j < TetPolyOrderDim[PolyFlag]; j++) { int index1 = tet_id * (TetPolyOrderDim[PolyFlag] * TetPolyOrderDim[PolyFlag]) + j * TetPolyOrderDim[PolyFlag] + i; int index2 = tet_id * TetPolyOrderDim[PolyFlag] + i; sum += InvMat[index1] * nxe_inc[j] * map[index2]; /* if (index2 == 0) { printf("INDEX1 = %d | INDEX2 = %d | map[index2] = %d | nxe_inc = %0.6e\n",index1,index2,map[index2],nxe_inc[j]); } */ } Field[tet_id * TetPolyOrderDim[PolyFlag] + i] += sum; /* if (ExcitationFacesOffset[tet_id] == 7) { printf("SUM = %0.6e | ind = %d \n",sum, tet_id * TetPolyOrderDim[PolyFlag] + i); } */ } } } __global__ void checker(fp_t_ts* a, fp_t_ts* b, int length){ int idx = blockIdx.x * blockDim.x + threadIdx.x; if(idx < length){ if(a[idx] != b[idx]){ printf("ERROR %.10f\n", a[idx] - b[idx]); } } } // ---- Small vector helpers (you probably already have similar ones) ---- __device__ inline void proj_to_plane(const fp_t_ts p[3], const fp_t_ts n[3], fp_t_ts out[3]) { // out = p - (p·n) n fp_t_ts dn = p[0]*n[0] + p[1]*n[1] + p[2]*n[2]; out[0] = p[0] - dn*n[0]; out[1] = p[1] - dn*n[1]; out[2] = p[2] - dn*n[2]; } __device__ inline fp_t_ts dot3(const fp_t_ts a[3], const fp_t_ts b[3]) { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; } __device__ inline void add3(const fp_t_ts a[3], const fp_t_ts b[3], fp_t_ts out[3]) { out[0]=a[0]+b[0]; out[1]=a[1]+b[1]; out[2]=a[2]+b[2]; } __device__ inline void sub3(const fp_t_ts a[3], const fp_t_ts b[3], fp_t_ts out[3]) { out[0]=a[0]-b[0]; out[1]=a[1]-b[1]; out[2]=a[2]-b[2]; } __device__ inline fp_t_ts norm3(const fp_t_ts a[3]) { return sqrtf(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]); } __device__ inline void normed(const fp_t_ts a[3], fp_t_ts out[3]) { fp_t_ts n = norm3(a); if (n > (fp_t_ts)0) { out[0]=a[0]/n; out[1]=a[1]/n; out[2]=a[2]/n; } else { out[0]=0; out[1]=0; out[2]=0; } } // -------------------------------------------------------------------------------------- // Antenna Port Update // Implement the port excitation using analytical solution // TEM Port // -------------------------------------------------------------------------------------- HD void TimeModulationInc(int TimeDistributionFlag, fp_t_ts t, fp_t_ts to, fp_t_ts freq_m, fp_t_ts Emag, fp_t_ts tau, fp_t_ts CHIRP_BW_MHZ, fp_t_ts &IncidExcit) { fp_t_ts Exponent; fp_t_ts omega = 2.0 * Pi * freq_m * MEGA; fp_t_ts Neuman; fp_t_ts SinModul; fp_t_ts kr; fp_t_ts work[3]; switch (TimeDistributionFlag){ // Time Harmonic case 0: IncidExcit = Emag * sin(- omega * (t)); break; // Gauss Pulse case 1: Exponent = t - to; SinModul = cos(omega * (t - to)); IncidExcit = Emag * SinModul * exp(-(Exponent * Exponent) / (tau * tau)); break; // Neuman Pulse case 2: Exponent = t - to; Neuman = (2.0 * Exponent) / (tau * tau); IncidExcit = (Emag * Neuman) * exp(-(Exponent * Exponent) / (tau * tau)); break; // 3) Linear FM chirp starting at 'to', duration 'tau', bandwidth CHIRP_BW_MHZ, // Hann-tapered on [0, tau] to reduce spectral sidelobes. case 3: { // Time since chirp start (seconds) fp_t_ts tt = t - to; // Gate to [0, tau] if (tt >= (fp_t_ts)0.0 && tt <= tau) { // Sweep parameters const fp_t_ts B = (fp_t_ts)(CHIRP_BW_MHZ * MEGA); // Hz const fp_t_ts fc = (fp_t_ts)(freq_m * MEGA); // Hz const fp_t_ts f0 = fc - (fp_t_ts)0.5 * B; // start freq [Hz] const fp_t_ts k = (tau != (fp_t_ts)0 ? B / tau : (fp_t_ts)0); // Hz/s // Phase: 2π( f0*tt + 0.5*k*tt^2 ) const fp_t_ts phase = (fp_t_ts)(2.0 * Pi) * (f0 * tt + (fp_t_ts)0.5 * k * tt * tt); // Hann window on [0, tau]: w = 0.5*(1 - cos(2π*tt/tau)) const fp_t_ts w = (fp_t_ts)0.5 * ((fp_t_ts)1.0 - cos((fp_t_ts)(2.0 * Pi) * (tt / tau))); const fp_t_ts s = sin(phase) * w; IncidExcit = Emag * s; } else { // Outside the chirp window IncidExcit = (fp_t_ts)0.0; } } break; // 4) Linear FM chirp starting at 'to', duration 'tau', bandwidth CHIRP_BW_MHZ, case 4: { // Time since chirp start (seconds) fp_t_ts tt = t - to; // Gate to [0, tau] if (tt >= (fp_t_ts)0.0 && tt <= tau) { // Sweep parameters const fp_t_ts B = (fp_t_ts)(CHIRP_BW_MHZ * MEGA); // Hz const fp_t_ts fc = (fp_t_ts)(freq_m * MEGA); // Hz const fp_t_ts f0 = fc - (fp_t_ts)0.5 * B; // start freq [Hz] const fp_t_ts k = (tau != (fp_t_ts)0 ? B / tau : (fp_t_ts)0); // Hz/s // Phase: 2π( f0*tt + 0.5*k*tt^2 ) const fp_t_ts phase = (fp_t_ts)(2.0 * Pi) * (f0 * tt + (fp_t_ts)0.5 * k * tt * tt); const fp_t_ts s = sin(phase); IncidExcit = Emag * s; } else { // Outside the chirp window IncidExcit = (fp_t_ts)0.0; } } break; default: break; } } // Computes nxH_inc and nxnxE_inc at a quadrature point. // Inputs: // t,to,tau,freq_m,TimeDistributionFlag,CHIRP_BW_MHZ : time envelope // exci_prop : geometry & parameters for the selected port type // n : port unit normal (propagation dir), length=1 // t1,t2 : orthonormal tangential basis on the port plane // r : 3D position of the current quadrature point // // Outputs: // nxHinc[3] = n × H_inc // nxnxEinc[3] = n × (n × E_inc) __device__ void updateinc_PORT(fp_t_ts t, const ExcitationProp& exci_prop, const fp_t_ts r[3], const fp_t_ts *normalvtr, fp_t_ts* nxField, fp_t_ts* nxnxField, int mode) { fp_t_ts t1[3] = {exci_prop.t1[0],exci_prop.t1[1],exci_prop.t1[2]}; fp_t_ts t2[3] = {exci_prop.t2[0],exci_prop.t2[1],exci_prop.t2[2]}; fp_t_ts nn[3] = {exci_prop.PortDirection[0],exci_prop.PortDirection[1],exci_prop.PortDirection[2]}; fp_t_ts n[3]; normed(nn, n); // ===== 0) Constants & medium ===== const fp_t_ts pi = (fp_t_ts)3.14159265358979323846; // ----- 1) Time envelope (scalar) ----- fp_t_ts IncidExcit = (fp_t_ts)1.0; fp_t_ts Emag = (fp_t_ts)exci_prop.Emagnitude; TimeModulationInc(exci_prop.TimeDistributionFlag, t, exci_prop.to, exci_prop.freq_m, Emag, exci_prop.tau, exci_prop.CHIRP_BW_MHZ, IncidExcit); // ----- 2) Build tangential E_inc shape per PortFlag ----- fp_t_ts E_tan[3] = {0,0,0}; // TE_mn geometry parameters (only meaningful for PortFlag==3) const int m = exci_prop.m; const int nmn = exci_prop.n; // avoid name clash with normal 'n' const fp_t_ts a_rect = (exci_prop.rect_a > (fp_t_ts)0 ? exci_prop.rect_a : (fp_t_ts)1e-12); const fp_t_ts b_rect = (exci_prop.rect_b > (fp_t_ts)0 ? exci_prop.rect_b : (fp_t_ts)1e-12); // Coax parameters (only meaningful for PortFlag==2) fp_t_ts r0[3] = {exci_prop.r0_port[0], exci_prop.r0_port[1], exci_prop.r0_port[2]}; fp_t_ts r1[3] = {exci_prop.r1_port[0], exci_prop.r1_port[1], exci_prop.r1_port[2]}; fp_t_ts r2[3] = {exci_prop.r2_port[0], exci_prop.r2_port[1], exci_prop.r2_port[2]}; if (exci_prop.PortFlag == 1) { // ================= TEM rectangular ================= // E is uniform along vpath projected into the port plane. fp_t_ts ehat[3]; normed(exci_prop.vpath, ehat); E_tan[0] = ehat[0]*Emag*IncidExcit; E_tan[1] = ehat[1]*Emag*IncidExcit; E_tan[2] = ehat[2]*Emag*IncidExcit; //printf("Etan = %0.6f,%0.6f,%0.6f\n", E_tan[0], E_tan[1], E_tan[2]); } else if (exci_prop.PortFlag == 2) { // ================= TEM coax ========================= // E(r) = (V / (rho * ln(b/a))) * rho_hat (in port plane) // Use Emag as V (scales overall level). Region: a <= rho <= b. // --- TEM coax: E = (V / (rho ln(b/a))) * rho_hat on a<=rho<=b --- fp_t_ts v10[3], v20[3]; sub3(r1,r0,v10); sub3(r2,r0,v20); fp_t_ts v10p[3], v20p[3]; proj_to_plane(v10,n,v10p); proj_to_plane(v20,n,v20p); const fp_t_ts a = norm3(v10p); const fp_t_ts b = norm3(v20p); fp_t_ts rp0[3], s[3]; sub3(r, r0, rp0); proj_to_plane(rp0, n, s); const fp_t_ts rho = norm3(s); if (rho >= a && rho <= b && a > (fp_t_ts)0 && b > a) { fp_t_ts sr_hat[3]; normed(s, sr_hat); const fp_t_ts amp = Emag*IncidExcit / (rho * logf(b/a)); E_tan[0] = sr_hat[0]*amp; E_tan[1] = sr_hat[1]*amp; E_tan[2] = sr_hat[2]*amp; } } else { // ================= TE_mn rectangular (PortFlag == 3) ============== // Local coords: u in [0,a] along t1, v in [0,b] along t2 from uv0. // Shape (time-domain real form, scaled by Emag*IncidExcit): // E_u ∝ (nπ/b) * cos(mπ u/a) * sin(nπ v/b) // E_v ∝ −(mπ/a) * sin(mπ u/a) * cos(nπ v/b) // (signs chosen to be consistent with ∇·E_t = 0; overall scale via Emag) // H_t for BC is supplied by Zport you set; here we’ll just use H = (1/Zport) n×E. // Local (u,v) from uv0 fp_t_ts uv0[3] = {exci_prop.uv0[0], exci_prop.uv0[1], exci_prop.uv0[2]}; fp_t_ts drv[3]; sub3(r, uv0, drv); const fp_t_ts u = dot3(drv, t1); const fp_t_ts v = dot3(drv, t2); const fp_t_ts cu = cosf((fp_t_ts)m * pi * u / a_rect); const fp_t_ts su = sinf((fp_t_ts)m * pi * u / a_rect); const fp_t_ts cv = cosf((fp_t_ts)nmn * pi * v / b_rect); const fp_t_ts sv = sinf((fp_t_ts)nmn * pi * v / b_rect); const fp_t_ts mu_u = (nmn!=0) ? ((fp_t_ts)nmn*pi/b_rect) : (fp_t_ts)0; const fp_t_ts mu_v = (m !=0) ? ((fp_t_ts)m *pi/a_rect) : (fp_t_ts)0; fp_t_ts Eu = ( mu_u ) * ( cu * sv ) * (Emag*IncidExcit); fp_t_ts Ev = -( mu_v ) * ( su * cv ) * (Emag*IncidExcit); // back to 3D E_tan[0] = Eu*t1[0] + Ev*t2[0]; E_tan[1] = Eu*t1[1] + Ev*t2[1]; E_tan[2] = Eu*t1[2] + Ev*t2[2]; } // ----- 3) Zport (implicit calc if exci_prop.PortImpedance <= 0) ----- fp_t_ts Zport = exci_prop.PortImpedance; //printf("Zport = %0.6e\n",Zport); // ===== 4) One cross, then reuse for both outputs ===== fp_t_ts H_tan[3]; // Calculate Magnetic Field crossP3(n, E_tan, H_tan); const fp_t_ts invZ = 1.0 / Zport; scale3(H_tan, invZ); fp_t_ts work[3]; if (mode == 1) { // nxHinc crossP3(normalvtr, H_tan, nxField); // nxnxEinc crossP3(normalvtr, E_tan, work); crossP3(normalvtr, work, nxnxField); } // nxnxHinc and nxEinc else { // nxEinc crossP3(normalvtr, E_tan, nxField); // nxnxHinc crossP3(normalvtr, H_tan, work); crossP3(normalvtr, work, nxnxField); } // cleanup for (int i=0;i<3;++i) { if (fabs(nxField[i]) < (fp_t_ts)1e-12) nxField[i] = 0.0; if (fabs(nxnxField[i]) < (fp_t_ts)1e-12) nxnxField[i] = 0.0; } } // Port-only version of H-field excitation __global__ void addExcitationH_port( int* ExcitationFacesCnt, int* ExcitationFacesOffset, int* ExcitationFacesNum, int excitationCNT, int8_t* map, ExcitationProp* exci_prop_vec, int* PortFacePidx_d, int PolyFlag, fp_t_ts dt_muo, fp_t_ts t, fp_t_ts* nd_coords_tet, fp_t_ts* nd_coords_face, fp_t_ts* InvMat, fp_t_ts* Field) { int tet_id = blockIdx.x * blockDim.x + threadIdx.x; if (tet_id < excitationCNT) { //printf("ExcitationFacesOffset = %d\n",ExcitationFacesOffset[tet_id]); const int dim = TetPolyOrderDim[PolyFlag]; const int fdim = FacePolyOrderDim[PolyFlag]; const int Q = GAUSS_POINT_NUM[PolyFlag]; //-----TODO: Make this part dynamic---------------------------- const int TetPolyOrderDim2 = TetPolyOrderDim[3]; const int FacePolyOrderDim2 = FacePolyOrderDim[3]; fp_t_ts nxe_inc[TetPolyOrderDim2] = {0}; fp_t_ts nxnxh_inc[TetPolyOrderDim2] = {0}; //------------------------------------------------------------- fp_t_ts Z = 0.0; for (int i = 0; i < ExcitationFacesCnt[tet_id]; i++) { //printf("i = %d | limit = %d \n",i,ExcitationFacesCnt[tet_id]); int face_id = ExcitationFacesOffset[tet_id] + i; //printf("face_id = %d\n",face_id); int port_id = PortFacePidx_d[face_id]; ExcitationProp exci_prop = exci_prop_vec[port_id]; Z = exci_prop.PortImpedance; fp_t_ts nvtr[2][3]; fp_t_ts normvtr[3] = {0}; fp_t_ts area = 0; fp_t_ts nxentry; fp_t_ts nxnxentry; fp_t_ts z[3], w[3], r[3]; fp_t_ts wgt; //-----------------TODO: Make this part dynamic-------------------------------------------- fp_t_ts nxelem[FacePolyOrderDim2] = {0}; fp_t_ts nxnxelem[FacePolyOrderDim2] = {0}; //----------------------------------------------------------------------------------------- fp_t_ts nxE_inc_[3]; fp_t_ts nxnxH_inc_[3]; /* if (face_id == 7) { printf("%0.6e, %0.6e, %0.6e %0.6e, %0.6e, %0.6e %0.6e, %0.6e, %0.6e\n", nd_coords_face[3 * 3 * face_id],nd_coords_face[3 * 3 * face_id+1],nd_coords_face[3 * 3 * face_id+2], nd_coords_face[3 * 3 * face_id+3],nd_coords_face[3 * 3 * face_id+4],nd_coords_face[3 * 3 * face_id+5], nd_coords_face[3 * 3 * face_id+6],nd_coords_face[3 * 3 * face_id+7],nd_coords_face[3 * 3 * face_id+8]); } */ geometry2(&nd_coords_face[3 * 3 * face_id], &nd_coords_tet[4 * 3 * tet_id], normvtr, nvtr, &area); for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++) { nxentry = 0; nxnxentry = 0; for (int k = 0; k < Q; k++) { quadObtention(Q, k, z[0], z[1], z[2], wgt); faceBasisW(&nd_coords_face[3 * 3 * face_id], z, j, w); getCartesianCoord(&nd_coords_face[3 * 3 * face_id], z, r); updateinc_PORT(t, exci_prop, r, normvtr, nxE_inc_, nxnxH_inc_, 2); if (!(isnan(nxE_inc_[0]))) { nxentry += area * wgt * dotP3131(w, nxE_inc_); nxnxentry += area * wgt * dotP3131(w, nxnxH_inc_); } /* if (tet_id == 6) { printf("area = %0.6e | wgt = %0.6e | w = %0.6e,%0.6e,%0.6e | nxE_inc = %0.6e,%0.6e,%0.6e | nxnxH_inc = %0.6e,%0.6e,%0.6e \n", area,wgt, w[0], w[1], w[2],nxE_inc_[0],nxE_inc_[1],nxE_inc_[2],nxnxH_inc_[0],nxnxH_inc_[1],nxnxH_inc_[2]); } */ } nxelem[j] += nxentry; nxnxelem[j] += nxnxentry; } for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++) { // printf("%d\n", ExcitationFacesNum[face_id]); nxe_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxelem[j]; nxnxh_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxnxelem[j]; } } for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++) { /* if (tet_id == 3) { printf("tet=%d i=%d | nxe_inc=%.6e | dt_muo=%.6e Z=%.6e | nxnxh=%.6e | map[%d]=%d\n", tet_id, i, (double)nxe_inc[i], (double)dt_muo, (double)Z, (double)nxnxh_inc[i], map[tet_id * TetPolyOrderDim[PolyFlag] + i]); } */ nxe_inc[i] *= -dt_muo; nxe_inc[i] += (-Z * dt_muo) * nxnxh_inc[i]; nxe_inc[i] *= map[tet_id * TetPolyOrderDim[PolyFlag] + i]; } fp_t_ts sum; for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++){ sum = 0.0; for (int j = 0; j < TetPolyOrderDim[PolyFlag]; j++) { int index1 = tet_id * (TetPolyOrderDim[PolyFlag] * TetPolyOrderDim[PolyFlag]) + j * TetPolyOrderDim[PolyFlag] + i; int index2 = tet_id * TetPolyOrderDim[PolyFlag] + i; sum += InvMat[index1] * nxe_inc[j] * map[index2]; /* if (index2 == 0) { printf("INDEX1 = %d | INDEX2 = %d | map[index2] = %d | nxe_inc = %0.6e\n",index1,index2,map[index2],nxe_inc[j]); } */ } Field[tet_id * TetPolyOrderDim[PolyFlag] + i] += sum; /* if (ExcitationFacesOffset[tet_id] == 7) { printf("SUM = %0.6e | ind = %d \n",sum, tet_id * TetPolyOrderDim[PolyFlag] + i); } */ } } } __global__ void addExcitationE_port( int* ExcitationFacesCnt, int* ExcitationFacesOffset, int* ExcitationFacesNum, int excitationCNT, int8_t* map, ExcitationProp* exci_prop_vec, int* PortFacePidx_d, int PolyFlag, fp_t_ts dt_epsilon, fp_t_ts t, fp_t_ts* nd_coords_tet, fp_t_ts* nd_coords_face, fp_t_ts* InvMat, fp_t_ts* Field ) { int tet_id = blockIdx.x * blockDim.x + threadIdx.x; if (tet_id < excitationCNT) { //-----TODO: Make this part dynamic---------------------------- const int TetPolyOrderDim2 = TetPolyOrderDim[3]; const int FacePolyOrderDim2 = FacePolyOrderDim[3]; fp_t_ts nxh_inc[TetPolyOrderDim2] = {0}; fp_t_ts nxnxe_inc[TetPolyOrderDim2] = {0}; //------------------------------------------------------------- fp_t_ts Y = 0.0; for (int i = 0; i < ExcitationFacesCnt[tet_id]; i++) { int face_id = ExcitationFacesOffset[tet_id] + i; int port_id = PortFacePidx_d[face_id]; ExcitationProp exci_prop = exci_prop_vec[port_id]; Y = 1.0 / exci_prop.PortImpedance; fp_t_ts nvtr[2][3]; fp_t_ts normvtr[3] = {0}; fp_t_ts area = 0; fp_t_ts nxentry; fp_t_ts nxnxentry; fp_t_ts z[3], w[3], r[3]; fp_t_ts wgt; //-----------------TODO: Make this part dynamic-------------------------------------------- fp_t_ts nxelem[FacePolyOrderDim2] = {0}; fp_t_ts nxnxelem[FacePolyOrderDim2] = {0}; //----------------------------------------------------------------------------------------- fp_t_ts nxH_inc_[3]; fp_t_ts nxnxE_inc_[3]; geometry2(&nd_coords_face[3 * 3 * face_id], &nd_coords_tet[4 * 3 * tet_id], normvtr, nvtr, &area); const int Q = GAUSS_POINT_NUM[PolyFlag]; for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++) { nxentry = 0; nxnxentry = 0; for (int k = 0; k < Q; k++) { quadObtention(Q, k, z[0], z[1], z[2], wgt); faceBasisW(&nd_coords_face[3 * 3 * face_id], z, j, w); getCartesianCoord(&nd_coords_face[3 * 3 * face_id], z, r); updateinc_PORT(t, exci_prop, r, normvtr, nxH_inc_, nxnxE_inc_, 1); if (!(isnan(nxH_inc_[0]))) { nxentry += area * wgt * dotP3131(w, nxH_inc_); nxnxentry += area * wgt * dotP3131(w, nxnxE_inc_); } } nxelem[j] += nxentry; nxnxelem[j] += nxnxentry; } for (int j = 0; j < FacePolyOrderDim[PolyFlag]; j++) { nxh_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxelem[j]; nxnxe_inc[FAC2TET[ExcitationFacesNum[face_id]][j]] += nxnxelem[j]; } } for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++) { nxh_inc[i] *= dt_epsilon; nxh_inc[i] += (-Y * dt_epsilon) * nxnxe_inc[i]; nxh_inc[i] *= map[tet_id * TetPolyOrderDim[PolyFlag] + i]; } fp_t_ts sum; for (int i = 0; i < TetPolyOrderDim[PolyFlag]; i++) { sum = 0.0; for (int j = 0; j < TetPolyOrderDim[PolyFlag]; j++) { sum += InvMat[tet_id * (TetPolyOrderDim[PolyFlag] * TetPolyOrderDim[PolyFlag]) + j * TetPolyOrderDim[PolyFlag] + i] * nxh_inc[j] * map[tet_id * TetPolyOrderDim[PolyFlag] + i]; } Field[tet_id * TetPolyOrderDim[PolyFlag] + i] += sum; //printf("Sum = %0.6e \n",sum); } } } #endif #endif // DGTD_KERNELS_CU