From 7488dbe4b6da44e8ead17e95b57025352d6bef8a Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sat, 27 Aug 2022 16:40:36 +0200 Subject: fill_eq_coeffs_template: reindent --- src/fill_eq_coeffs_template.c | 934 +++++++++++++++++++++--------------------- 1 file changed, 467 insertions(+), 467 deletions(-) diff --git a/src/fill_eq_coeffs_template.c b/src/fill_eq_coeffs_template.c index ae8fcd8..20706f6 100644 --- a/src/fill_eq_coeffs_template.c +++ b/src/fill_eq_coeffs_template.c @@ -110,26 +110,26 @@ static void fill_eq_coeffs_line(const cGH *gh, const CoordPatch *cp, MG2DContext const double * const vars[], const int idx_z, const ptrdiff_t stride_z, const double dx_inv, const double dz_inv) { - for (int idx_x = 0; idx_x < solver->local_size[0]; idx_x++) { - const int idx_src = CCTK_GFINDEX3D(gh, idx_x + cp->offset_left[0], cp->y_idx, idx_z + cp->offset_left[1]); - const int idx_dc = idx_z * solver->diff_coeffs[0]->stride + idx_x; - const int idx_rhs = idx_z * solver->rhs_stride + idx_x; - - const double x = vars[RHS_VAR_X][idx_src]; - const int on_axis = fabs(x) < 1e-13; - - double K[3][3], Km[3][3], Ku[3][3], dK[3][3][3]; - double gtu[3][3], g[3][3], gu[3][3]; - double dg[3][3][3], dgu[3][3][3], d2g[3][3][3][3], gudot[3][3]; - double G[3][3][3], dG[3][3][3][3], Gdot[3][3][3], X[3]; - double Ric[3][3], Ricm[3][3]; - double dgdot[3][3][3]; - double D2alpha[3][3]; - double Kdot[3][3]; - double k2, Kij_Dij_alpha, k_kdot, k3; - double Gammadot_term, beta_term; - - double betal[3], dbetal[3][3], d2betal[3][3][3]; + for (int idx_x = 0; idx_x < solver->local_size[0]; idx_x++) { + const int idx_src = CCTK_GFINDEX3D(gh, idx_x + cp->offset_left[0], cp->y_idx, idx_z + cp->offset_left[1]); + const int idx_dc = idx_z * solver->diff_coeffs[0]->stride + idx_x; + const int idx_rhs = idx_z * solver->rhs_stride + idx_x; + + const double x = vars[RHS_VAR_X][idx_src]; + const int on_axis = fabs(x) < 1e-13; + + double K[3][3], Km[3][3], Ku[3][3], dK[3][3][3]; + double gtu[3][3], g[3][3], gu[3][3]; + double dg[3][3][3], dgu[3][3][3], d2g[3][3][3][3], gudot[3][3]; + double G[3][3][3], dG[3][3][3][3], Gdot[3][3][3], X[3]; + double Ric[3][3], Ricm[3][3]; + double dgdot[3][3][3]; + double D2alpha[3][3]; + double Kdot[3][3]; + double k2, Kij_Dij_alpha, k_kdot, k3; + double Gammadot_term, beta_term; + + double betal[3], dbetal[3][3], d2betal[3][3][3]; #define LOAD_VAR(var) ( vars[RHS_VAR_ ## var][idx_src]) #define DX(var) (diff1 (&vars[RHS_VAR_ ## var][idx_src], 1, dx_inv )) @@ -138,217 +138,217 @@ static void fill_eq_coeffs_line(const cGH *gh, const CoordPatch *cp, MG2DContext #define DZZ(var) (diff2 (&vars[RHS_VAR_ ## var][idx_src], stride_z, dz_inv)) #define DXZ(var) (diff11(&vars[RHS_VAR_ ## var][idx_src], 1, stride_z, dx_inv, dz_inv)) - const double gtxx = LOAD_VAR(GTXX); - const double gtyy = LOAD_VAR(GTYY); - const double gtzz = LOAD_VAR(GTZZ); - const double gtxy = LOAD_VAR(GTXY); - const double gtxz = LOAD_VAR(GTXZ); - const double gtyz = LOAD_VAR(GTYZ); - - const double gt[3][3] = {{ gtxx, gtxy, gtxz }, - { gtxy, gtyy, gtyz }, - { gtxz, gtyz, gtzz }}; - - const double dx_gt11 = DX(GTXX); - const double dx_gt22 = DX(GTYY); - const double dx_gt33 = DX(GTZZ); - const double dx_gt13 = DX(GTXZ); - - const double dz_gt11 = DZ(GTXX); - const double dz_gt22 = DZ(GTYY); - const double dz_gt33 = DZ(GTZZ); - const double dz_gt13 = DZ(GTXZ); - - const double dgt[3][3][3] = { + const double gtxx = LOAD_VAR(GTXX); + const double gtyy = LOAD_VAR(GTYY); + const double gtzz = LOAD_VAR(GTZZ); + const double gtxy = LOAD_VAR(GTXY); + const double gtxz = LOAD_VAR(GTXZ); + const double gtyz = LOAD_VAR(GTYZ); + + const double gt[3][3] = {{ gtxx, gtxy, gtxz }, + { gtxy, gtyy, gtyz }, + { gtxz, gtyz, gtzz }}; + + const double dx_gt11 = DX(GTXX); + const double dx_gt22 = DX(GTYY); + const double dx_gt33 = DX(GTZZ); + const double dx_gt13 = DX(GTXZ); + + const double dz_gt11 = DZ(GTXX); + const double dz_gt22 = DZ(GTYY); + const double dz_gt33 = DZ(GTZZ); + const double dz_gt13 = DZ(GTXZ); + + const double dgt[3][3][3] = { + { + { dx_gt11, 0.0, dx_gt13 }, + { 0.0, dx_gt22, 0.0 }, + { dx_gt13, 0.0, dx_gt33 }, + }, + { + { 0.0, on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0 }, + { on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0, on_axis ? dx_gt13 : gtxz / x }, + { 0.0, on_axis ? dx_gt13 : gtxz / x, 0.0 }, + }, + { + { dz_gt11, 0.0, dz_gt13 }, + { 0.0, dz_gt22, 0.0 }, + { dz_gt13, 0.0, dz_gt33 }, + }, + }; + + const double dxx_gt11 = DXX(GTXX); + const double dxx_gt22 = DXX(GTYY); + const double dxx_gt33 = DXX(GTZZ); + const double dxx_gt13 = DXX(GTXZ); + + const double dxz_gt11 = DXZ(GTXX); + const double dxz_gt22 = DXZ(GTYY); + const double dxz_gt33 = DXZ(GTZZ); + const double dxz_gt13 = DXZ(GTXZ); + + const double dzz_gt11 = DZZ(GTXX); + const double dzz_gt22 = DZZ(GTYY); + const double dzz_gt33 = DZZ(GTZZ); + const double dzz_gt13 = DZZ(GTXZ); + + const double d2gt[3][3][3][3] = { + { { - { dx_gt11, 0.0, dx_gt13 }, - { 0.0, dx_gt22, 0.0 }, - { dx_gt13, 0.0, dx_gt33 }, + { dxx_gt11, 0.0, dxx_gt13 }, + { 0.0, dxx_gt22, 0.0 }, + { dxx_gt13, 0.0, dxx_gt33 }, }, { - { 0.0, on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0 }, - { on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0, on_axis ? dx_gt13 : gtxz / x }, - { 0.0, on_axis ? dx_gt13 : gtxz / x, 0.0 }, + { 0.0, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 }, + { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0, + on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, + { 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 }, }, { - { dz_gt11, 0.0, dz_gt13 }, - { 0.0, dz_gt22, 0.0 }, - { dz_gt13, 0.0, dz_gt33 }, + { dxz_gt11, 0.0, dxz_gt13 }, + { 0.0, dxz_gt22, 0.0 }, + { dxz_gt13, 0.0, dxz_gt33 }, }, - }; - const double dxx_gt11 = DXX(GTXX); - const double dxx_gt22 = DXX(GTYY); - const double dxx_gt33 = DXX(GTZZ); - const double dxx_gt13 = DXX(GTXZ); - - const double dxz_gt11 = DXZ(GTXX); - const double dxz_gt22 = DXZ(GTYY); - const double dxz_gt33 = DXZ(GTZZ); - const double dxz_gt13 = DXZ(GTXZ); - - const double dzz_gt11 = DZZ(GTXX); - const double dzz_gt22 = DZZ(GTYY); - const double dzz_gt33 = DZZ(GTZZ); - const double dzz_gt13 = DZZ(GTXZ); - - const double d2gt[3][3][3][3] = { + }, + { { - { - { dxx_gt11, 0.0, dxx_gt13 }, - { 0.0, dxx_gt22, 0.0 }, - { dxx_gt13, 0.0, dxx_gt33 }, - }, - { - { 0.0, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 }, - { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0, - on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, - { 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 }, - }, - { - { dxz_gt11, 0.0, dxz_gt13 }, - { 0.0, dxz_gt22, 0.0 }, - { dxz_gt13, 0.0, dxz_gt33 }, - }, - + { 0.0, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 }, + { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0, + on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, + { 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 }, }, { - { - { 0.0, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 }, - { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0, - on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, - { 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 }, - }, - { - { on_axis ? dxx_gt22 : dx_gt11 / x - 2 * (gtxx - gtyy) / SQR(x), 0.0, - on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, - { 0.0, on_axis ? dxx_gt11 : dx_gt22 / x + 2.0 * (gtxx - gtyy) / SQR(x), 0.0 }, - { on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0, on_axis ? dxx_gt33 : dx_gt33 / x }, - }, - { - { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 }, - { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0, - on_axis ? dxz_gt13 : dz_gt13 / x }, - { 0.0, on_axis ? dxz_gt13 : dz_gt13 / x, 0.0 }, - }, - + { on_axis ? dxx_gt22 : dx_gt11 / x - 2 * (gtxx - gtyy) / SQR(x), 0.0, + on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, + { 0.0, on_axis ? dxx_gt11 : dx_gt22 / x + 2.0 * (gtxx - gtyy) / SQR(x), 0.0 }, + { on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0, on_axis ? dxx_gt33 : dx_gt33 / x }, }, { - { - { dxz_gt11, 0.0, dxz_gt13 }, - { 0.0, dxz_gt22, 0.0 }, - { dxz_gt13, 0.0, dxz_gt33 }, - }, - { - { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 }, - { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0, - on_axis ? dxz_gt13 : dz_gt13 / x }, - { 0.0, on_axis ? dxz_gt13 : dz_gt13 / x, 0.0 }, - }, - { - { dzz_gt11, 0.0, dzz_gt13 }, - { 0.0, dzz_gt22, 0.0 }, - { dzz_gt13, 0.0, dzz_gt33 }, - }, - + { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 }, + { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0, + on_axis ? dxz_gt13 : dz_gt13 / x }, + { 0.0, on_axis ? dxz_gt13 : dz_gt13 / x, 0.0 }, }, - }; - - const double Atxx = LOAD_VAR(ATXX); - const double Atyy = LOAD_VAR(ATYY); - const double Atzz = LOAD_VAR(ATZZ); - const double Atxy = LOAD_VAR(ATXY); - const double Atxz = LOAD_VAR(ATXZ); - const double Atyz = LOAD_VAR(ATYZ); - - const double At[3][3] = { - { Atxx, Atxy, Atxz }, - { Atxy, Atyy, Atyz }, - { Atxz, Atyz, Atzz }}; - - const double dx_At11 = DX(ATXX); - const double dx_At22 = DX(ATYY); - const double dx_At33 = DX(ATZZ); - const double dx_At13 = DX(ATXZ); - - const double dz_At11 = DZ(ATXX); - const double dz_At22 = DZ(ATYY); - const double dz_At33 = DZ(ATZZ); - const double dz_At13 = DZ(ATXZ); - - const double dAt[3][3][3] = { + + }, + { { - { dx_At11, 0.0, dx_At13 }, - { 0.0, dx_At22, 0.0 }, - { dx_At13, 0.0, dx_At33 }, + { dxz_gt11, 0.0, dxz_gt13 }, + { 0.0, dxz_gt22, 0.0 }, + { dxz_gt13, 0.0, dxz_gt33 }, }, { - { 0.0, on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0 }, - { on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0, on_axis ? dx_At13 : Atxz / x }, - { 0.0, on_axis ? dx_At13 : Atxz / x, 0.0 }, + { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 }, + { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0, + on_axis ? dxz_gt13 : dz_gt13 / x }, + { 0.0, on_axis ? dxz_gt13 : dz_gt13 / x, 0.0 }, }, { - { dz_At11, 0.0, dz_At13 }, - { 0.0, dz_At22, 0.0 }, - { dz_At13, 0.0, dz_At33 }, + { dzz_gt11, 0.0, dzz_gt13 }, + { 0.0, dzz_gt22, 0.0 }, + { dzz_gt13, 0.0, dzz_gt33 }, }, - }; - - const double phi = LOAD_VAR(PHI); - const double dphi[3] = { DX(PHI), 0.0, DZ(PHI) }; - - const double dxx_phi = DXX(PHI); - const double dzz_phi = DZZ(PHI); - const double dxz_phi = DXZ(PHI); - - const double d2phi[3][3] = { - { dxx_phi, 0.0, dxz_phi }, - { 0.0, on_axis ? dxx_phi : dphi[0] / x, 0.0 }, - { dxz_phi, 0.0, dzz_phi }, - }; - - const double trK = LOAD_VAR(TRK); - const double dtrK[3] = { DX(TRK), 0.0, DZ(TRK) }; - const double Xtx = LOAD_VAR(XTX); - const double Xtz = LOAD_VAR(XTZ); - - const double alpha = LOAD_VAR(ALPHA); - const double dx_alpha = DX(ALPHA); - const double dz_alpha = DZ(ALPHA); - - const double dalpha[3] = { dx_alpha, 0.0, dz_alpha }; - - const double dxx_alpha = DXX(ALPHA); - const double dzz_alpha = DZZ(ALPHA); - const double dxz_alpha = DXZ(ALPHA); - - const double d2alpha[3][3] = {{ dxx_alpha, 0, dxz_alpha }, - { 0, on_axis ? dxx_alpha : dx_alpha / x, 0 }, - { dxz_alpha, 0, dzz_alpha }}; - - const double betax = LOAD_VAR(BETAX); - const double betaz = LOAD_VAR(BETAX); - - const double beta[3] = { betax, 0.0, betaz }; - - const double dx_betax = DX(BETAX); - const double dz_betax = DZ(BETAX); - - const double dx_betaz = DX(BETAZ); - const double dz_betaz = DZ(BETAZ); - - const double dbeta[3][3] = {{ dx_betax, 0.0, dx_betaz }, - { 0.0, on_axis ? dx_betax : betax / x, 0.0 }, - { dz_betax, 0.0, dz_betaz }}; - - const double dxx_betax = DXX(BETAX); - const double dxz_betax = DXZ(BETAX); - const double dzz_betax = DZZ(BETAX); - - const double dxx_betaz = DXX(BETAZ); - const double dxz_betaz = DXZ(BETAZ); - const double dzz_betaz = DZZ(BETAZ); + }, + }; + + const double Atxx = LOAD_VAR(ATXX); + const double Atyy = LOAD_VAR(ATYY); + const double Atzz = LOAD_VAR(ATZZ); + const double Atxy = LOAD_VAR(ATXY); + const double Atxz = LOAD_VAR(ATXZ); + const double Atyz = LOAD_VAR(ATYZ); + + const double At[3][3] = { + { Atxx, Atxy, Atxz }, + { Atxy, Atyy, Atyz }, + { Atxz, Atyz, Atzz }}; + + const double dx_At11 = DX(ATXX); + const double dx_At22 = DX(ATYY); + const double dx_At33 = DX(ATZZ); + const double dx_At13 = DX(ATXZ); + + const double dz_At11 = DZ(ATXX); + const double dz_At22 = DZ(ATYY); + const double dz_At33 = DZ(ATZZ); + const double dz_At13 = DZ(ATXZ); + + const double dAt[3][3][3] = { + { + { dx_At11, 0.0, dx_At13 }, + { 0.0, dx_At22, 0.0 }, + { dx_At13, 0.0, dx_At33 }, + }, + { + { 0.0, on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0 }, + { on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0, on_axis ? dx_At13 : Atxz / x }, + { 0.0, on_axis ? dx_At13 : Atxz / x, 0.0 }, + }, + { + { dz_At11, 0.0, dz_At13 }, + { 0.0, dz_At22, 0.0 }, + { dz_At13, 0.0, dz_At33 }, + }, + }; + + const double phi = LOAD_VAR(PHI); + const double dphi[3] = { DX(PHI), 0.0, DZ(PHI) }; + + const double dxx_phi = DXX(PHI); + const double dzz_phi = DZZ(PHI); + const double dxz_phi = DXZ(PHI); + + const double d2phi[3][3] = { + { dxx_phi, 0.0, dxz_phi }, + { 0.0, on_axis ? dxx_phi : dphi[0] / x, 0.0 }, + { dxz_phi, 0.0, dzz_phi }, + }; + + const double trK = LOAD_VAR(TRK); + const double dtrK[3] = { DX(TRK), 0.0, DZ(TRK) }; + + const double Xtx = LOAD_VAR(XTX); + const double Xtz = LOAD_VAR(XTZ); + + const double alpha = LOAD_VAR(ALPHA); + const double dx_alpha = DX(ALPHA); + const double dz_alpha = DZ(ALPHA); + + const double dalpha[3] = { dx_alpha, 0.0, dz_alpha }; + + const double dxx_alpha = DXX(ALPHA); + const double dzz_alpha = DZZ(ALPHA); + const double dxz_alpha = DXZ(ALPHA); + + const double d2alpha[3][3] = {{ dxx_alpha, 0, dxz_alpha }, + { 0, on_axis ? dxx_alpha : dx_alpha / x, 0 }, + { dxz_alpha, 0, dzz_alpha }}; + + const double betax = LOAD_VAR(BETAX); + const double betaz = LOAD_VAR(BETAX); + + const double beta[3] = { betax, 0.0, betaz }; + + const double dx_betax = DX(BETAX); + const double dz_betax = DZ(BETAX); + + const double dx_betaz = DX(BETAZ); + const double dz_betaz = DZ(BETAZ); + + const double dbeta[3][3] = {{ dx_betax, 0.0, dx_betaz }, + { 0.0, on_axis ? dx_betax : betax / x, 0.0 }, + { dz_betax, 0.0, dz_betaz }}; + + const double dxx_betax = DXX(BETAX); + const double dxz_betax = DXZ(BETAX); + const double dzz_betax = DZZ(BETAX); + + const double dxx_betaz = DXX(BETAZ); + const double dxz_betaz = DXZ(BETAZ); + const double dzz_betaz = DZZ(BETAZ); #undef LOAD_VAR #undef DX @@ -357,308 +357,308 @@ static void fill_eq_coeffs_line(const cGH *gh, const CoordPatch *cp, MG2DContext #undef DZZ #undef DXZ - const double d2beta[3][3][3] = { - { - { dxx_betax, 0.0, dxx_betaz }, - { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 }, - { dxz_betax, 0.0, dxz_betaz }, - }, - { - { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 }, - { on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0, on_axis ? dxx_betaz : dx_betaz / x }, - { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 }, - }, - { - { dxz_betax, 0.0, dxz_betaz }, - { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 }, - { dzz_betax, 0.0, dzz_betaz }, - }, - }; - - const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz); - - // \tilde{γ}^{ij} - gtu[0][0] = (gtyy * gtzz - SQR(gtyz)) / det; - gtu[1][1] = (gtxx * gtzz - SQR(gtxz)) / det; - gtu[2][2] = (gtxx * gtyy - SQR(gtxy)) / det; - gtu[0][1] = -(gtxy * gtzz - gtyz * gtxz) / det; - gtu[0][2] = (gtxy * gtyz - gtyy * gtxz) / det; - gtu[1][2] = -(gtxx * gtyz - gtxy * gtxz) / det; - gtu[1][0] = gtu[0][1]; - gtu[2][0] = gtu[0][2]; - gtu[2][1] = gtu[1][2]; - - //　γ_{jk}/^{jk} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { - gu[j][k] = SQR(phi) * gtu[j][k]; - g[j][k] = gt[j][k] / SQR(phi); - } - - // β_j - for (int j = 0; j < 3; j++) { - double val = 0.0; - for (int k = 0; k < 3; k++) - val += g[j][k] * beta[k]; - betal[j] = val; + const double d2beta[3][3][3] = { + { + { dxx_betax, 0.0, dxx_betaz }, + { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 }, + { dxz_betax, 0.0, dxz_betaz }, + }, + { + { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 }, + { on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0, on_axis ? dxx_betaz : dx_betaz / x }, + { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 }, + }, + { + { dxz_betax, 0.0, dxz_betaz }, + { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 }, + { dzz_betax, 0.0, dzz_betaz }, + }, + }; + + const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz); + + // \tilde{γ}^{ij} + gtu[0][0] = (gtyy * gtzz - SQR(gtyz)) / det; + gtu[1][1] = (gtxx * gtzz - SQR(gtxz)) / det; + gtu[2][2] = (gtxx * gtyy - SQR(gtxy)) / det; + gtu[0][1] = -(gtxy * gtzz - gtyz * gtxz) / det; + gtu[0][2] = (gtxy * gtyz - gtyy * gtxz) / det; + gtu[1][2] = -(gtxx * gtyz - gtxy * gtxz) / det; + gtu[1][0] = gtu[0][1]; + gtu[2][0] = gtu[0][2]; + gtu[2][1] = gtu[1][2]; + + //　γ_{jk}/^{jk} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + gu[j][k] = SQR(phi) * gtu[j][k]; + g[j][k] = gt[j][k] / SQR(phi); } - // ∂_j γ_{kl} - // ∂_j K_{kl} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) { - dg[j][k][l] = -2.0 * dphi[j] * gt[k][l] / (phi * SQR(phi)) + dgt[j][k][l] / SQR(phi); - dK[j][k][l] = -2.0 * dphi[j] * At[k][l] / (phi * SQR(phi)) + dAt[j][k][l] / SQR(phi) + - (1.0 / 3.0) * (dtrK[j] * g[k][l] + trK * dg[j][k][l]); - } - - // ∂_{jk} g_{lm} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) - for (int m = 0; m < 3; m++) { - d2g[j][k][l][m] = 6.0 * gt [l][m] * dphi[j] * dphi[k] / SQR(SQR(phi)) - - 2.0 * gt [l][m] * d2phi[j][k] / (phi * SQR(phi)) - - 2.0 * dgt [j][l][m] * dphi[k] / (phi * SQR(phi)) - - 2.0 * dgt [k][l][m] * dphi[j] / (phi * SQR(phi)) + - d2gt[j][k][l][m] / SQR(phi); - } + // β_j + for (int j = 0; j < 3; j++) { + double val = 0.0; + for (int k = 0; k < 3; k++) + val += g[j][k] * beta[k]; + betal[j] = val; + } - // ∂_j γ^{kl} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) { - double val = 0.0; - for (int m = 0; m < 3; m++) - for (int n = 0; n < 3; n++) - val += -gu[k][m] * gu[l][n] * dg[j][m][n]; - dgu[j][k][l] = val; - } + // ∂_j γ_{kl} + // ∂_j K_{kl} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) { + dg[j][k][l] = -2.0 * dphi[j] * gt[k][l] / (phi * SQR(phi)) + dgt[j][k][l] / SQR(phi); + dK[j][k][l] = -2.0 * dphi[j] * At[k][l] / (phi * SQR(phi)) + dAt[j][k][l] / SQR(phi) + + (1.0 / 3.0) * (dtrK[j] * g[k][l] + trK * dg[j][k][l]); + } - // Γ^j_{kl} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) { - double val = 0.0; - for (int m = 0; m < 3; m++) - val += 0.5 * gu[j][m] * (dg[k][l][m] + dg[l][k][m] - dg[m][k][l]); - G[j][k][l] = val; + // ∂_{jk} g_{lm} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + for (int m = 0; m < 3; m++) { + d2g[j][k][l][m] = 6.0 * gt [l][m] * dphi[j] * dphi[k] / SQR(SQR(phi)) - + 2.0 * gt [l][m] * d2phi[j][k] / (phi * SQR(phi)) - + 2.0 * dgt [j][l][m] * dphi[k] / (phi * SQR(phi)) - + 2.0 * dgt [k][l][m] * dphi[j] / (phi * SQR(phi)) + + d2gt[j][k][l][m] / SQR(phi); } - // ∂_j Γ^k_{lm} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) - for (int m = 0; m < 3; m++) { - double val = 0.0; - for (int n = 0; n < 3; n++) { - val += dgu[j][k][n] * (dg [l][m][n] + dg [m][l][n] - dg [n][l][m]) + - gu [k][n] * (d2g[j][l][m][n] + d2g[j][m][l][n] - d2g[j][n][l][m]); - } - dG[j][k][l][m] = 0.5 * val; - } - - // Γ^j = γ^{kl} Γ_{kl}^j - for (int j = 0; j < 3; j++) { - double val = 0.0; - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) - val += gu[k][l] * G[j][k][l]; - X[j] = val; - } + // ∂_j γ^{kl} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) { + double val = 0.0; + for (int m = 0; m < 3; m++) + for (int n = 0; n < 3; n++) + val += -gu[k][m] * gu[l][n] * dg[j][m][n]; + dgu[j][k][l] = val; + } - // ∂_j β_k - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { + // Γ^j_{kl} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) { double val = 0.0; - for (int l = 0; l < 3; l++) - val += g[k][l] * dbeta[j][l] + dg[j][k][l] * beta[l]; - dbetal[j][k] = val; + for (int m = 0; m < 3; m++) + val += 0.5 * gu[j][m] * (dg[k][l][m] + dg[l][k][m] - dg[m][k][l]); + G[j][k][l] = val; } - // ∂_{jk} β_l - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) { + // ∂_j Γ^k_{lm} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + for (int m = 0; m < 3; m++) { double val = 0.0; - for (int m = 0; m < 3; m++) - val += d2g[j][k][l][m] * beta[m] + dg[k][l][m] * dbeta[j][m] + dg[j][l][m] * dbeta[k][m] + - g[l][m] * d2beta[j][k][m]; - d2betal[j][k][l] = val; + for (int n = 0; n < 3; n++) { + val += dgu[j][k][n] * (dg [l][m][n] + dg [m][l][n] - dg [n][l][m]) + + gu [k][n] * (d2g[j][l][m][n] + d2g[j][m][l][n] - d2g[j][n][l][m]); + } + dG[j][k][l][m] = 0.5 * val; } - // K_{ij} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - K[j][k] = At[j][k] / SQR(phi) + (1.0 / 3.0) * g[j][k] * trK; + // Γ^j = γ^{kl} Γ_{kl}^j + for (int j = 0; j < 3; j++) { + double val = 0.0; + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + val += gu[k][l] * G[j][k][l]; + X[j] = val; + } - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { - double val = 0.0; - for (int l = 0; l < 3; l++) - val += gu[j][l] * K[l][k]; - Km[j][k] = val; - } + // ∂_j β_k + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + double val = 0.0; + for (int l = 0; l < 3; l++) + val += g[k][l] * dbeta[j][l] + dg[j][k][l] * beta[l]; + dbetal[j][k] = val; + } - // K^{ij} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { + // ∂_{jk} β_l + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) { double val = 0.0; - for (int l = 0; l < 3; l++) - val += gu[j][l] * Km[k][l]; - Ku[j][k] = val; + for (int m = 0; m < 3; m++) + val += d2g[j][k][l][m] * beta[m] + dg[k][l][m] * dbeta[j][m] + dg[j][l][m] * dbeta[k][m] + + g[l][m] * d2beta[j][k][m]; + d2betal[j][k][l] = val; } - // ∂_j \dot{γ_kl} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) { - dgdot[j][k][l] = -2.0 * (dalpha[j] * K[k][l] + alpha * dK[j][k][l]) + - d2betal[j][k][l] + d2betal[j][l][k]; - for (int m = 0; m < 3; m++) - dgdot[j][k][l] += -2.0 * (dG[j][m][k][l] * betal[m] + G[m][k][l] * dbetal[j][m]); - } + // K_{ij} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + K[j][k] = At[j][k] / SQR(phi) + (1.0 / 3.0) * g[j][k] * trK; - // \dot{γ^{jk}} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { - gudot[j][k] = 2.0 * alpha * Ku[j][k]; + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + double val = 0.0; + for (int l = 0; l < 3; l++) + val += gu[j][l] * K[l][k]; + Km[j][k] = val; + } - for (int l = 0; l < 3; l++) { - gudot[j][k] += -gu[j][l] * dbeta[l][k] - gu[k][l] * dbeta[l][j]; + // K^{ij} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + double val = 0.0; + for (int l = 0; l < 3; l++) + val += gu[j][l] * Km[k][l]; + Ku[j][k] = val; + } - for (int m = 0; m < 3; m++) - gudot[j][k] += -gu[j][l] * G[k][l][m] * beta[m] - gu[k][l] * G[j][l][m] * beta[m]; - } + // ∂_j \dot{γ_kl} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) { + dgdot[j][k][l] = -2.0 * (dalpha[j] * K[k][l] + alpha * dK[j][k][l]) + + d2betal[j][k][l] + d2betal[j][l][k]; + for (int m = 0; m < 3; m++) + dgdot[j][k][l] += -2.0 * (dG[j][m][k][l] * betal[m] + G[m][k][l] * dbetal[j][m]); } + // \dot{γ^{jk}} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + gudot[j][k] = 2.0 * alpha * Ku[j][k]; - // \dot{Γ}^j_{kl} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) { - double val = 0.0; - for (int m = 0; m < 3; m++) - val += gudot[j][m] * (dg [k][l][m] + dg [l][k][m] - dg [m][k][l]) + - gu [j][m] * (dgdot[k][l][m] + dgdot[l][k][m] - dgdot[m][k][l]); + for (int l = 0; l < 3; l++) { + gudot[j][k] += -gu[j][l] * dbeta[l][k] - gu[k][l] * dbeta[l][j]; - Gdot[j][k][l] = 0.5 * val; - } - - Gammadot_term = 0.0; - for (int j = 0; j < 3; j++) { - double val = 0.0; - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) { - val += gu[k][l] * Gdot[j][k][l]; - } - - Gammadot_term += val * dalpha[j]; + for (int m = 0; m < 3; m++) + gudot[j][k] += -gu[j][l] * G[k][l][m] * beta[m] - gu[k][l] * G[j][l][m] * beta[m]; + } } - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { + + // \dot{Γ}^j_{kl} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) { double val = 0.0; - for (int l = 0; l < 3; l++) - val += G[l][j][k] * dalpha[l]; + for (int m = 0; m < 3; m++) + val += gudot[j][m] * (dg [k][l][m] + dg [l][k][m] - dg [m][k][l]) + + gu [j][m] * (dgdot[k][l][m] + dgdot[l][k][m] - dgdot[m][k][l]); - D2alpha[j][k] = d2alpha[j][k] - val; + Gdot[j][k][l] = 0.5 * val; } - Kij_Dij_alpha = 0.0; - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { - Kij_Dij_alpha += Ku[j][k] * D2alpha[j][k]; + Gammadot_term = 0.0; + for (int j = 0; j < 3; j++) { + double val = 0.0; + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) { + val += gu[k][l] * Gdot[j][k][l]; } - k3 = 0.0; - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { - double val = 0.0; - for (int l = 0; l < 3; l++) - val += Km[k][l] * Km[l][j]; - k3 += val * Km[j][k]; - } + Gammadot_term += val * dalpha[j]; + } + + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + double val = 0.0; + for (int l = 0; l < 3; l++) + val += G[l][j][k] * dalpha[l]; - // K_{ij} K^{ij} - k2 = 0.0; - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - k2 += Km[j][k] * Km[k][j]; + D2alpha[j][k] = d2alpha[j][k] - val; + } - // Ric_{jk} - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { - double val = 0.0; - for (int m = 0; m < 3; m++) - val += dG[m][m][j][k] - dG[k][m][j][m]; - for (int m = 0; m < 3; m++) - for (int l = 0; l < 3; l++) - val += G[l][l][m] * G[m][j][k] - G[l][k][m] * G[m][j][l]; - Ric[j][k] = val; - } + Kij_Dij_alpha = 0.0; + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + Kij_Dij_alpha += Ku[j][k] * D2alpha[j][k]; + } - // Ric^j_k - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { - double val = 0.0; - for (int l = 0; l < 3; l++) - val += gu[j][l] * Ric[l][k]; - Ricm[j][k] = val; - } + k3 = 0.0; + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + double val = 0.0; + for (int l = 0; l < 3; l++) + val += Km[k][l] * Km[l][j]; + k3 += val * Km[j][k]; + } - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { - double val = 0.0; + // K_{ij} K^{ij} + k2 = 0.0; + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + k2 += Km[j][k] * Km[k][j]; + + // Ric_{jk} + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + double val = 0.0; + for (int m = 0; m < 3; m++) + val += dG[m][m][j][k] - dG[k][m][j][m]; + for (int m = 0; m < 3; m++) for (int l = 0; l < 3; l++) - val += K[j][l] * Km[l][k]; + val += G[l][l][m] * G[m][j][k] - G[l][k][m] * G[m][j][l]; + Ric[j][k] = val; + } - Kdot[j][k] = -D2alpha[j][k] + alpha * (Ric[j][k] + trK * K[j][k] - 2 * val); - } + // Ric^j_k + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + double val = 0.0; + for (int l = 0; l < 3; l++) + val += gu[j][l] * Ric[l][k]; + Ricm[j][k] = val; + } - k_kdot = 0.0; - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) { - //k_kdot += kdot[j][k] * Ku[j][k]; - k_kdot += Kdot[j][k] * Ku[j][k]; - } + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + double val = 0.0; + for (int l = 0; l < 3; l++) + val += K[j][l] * Km[l][k]; + + Kdot[j][k] = -D2alpha[j][k] + alpha * (Ric[j][k] + trK * K[j][k] - 2 * val); + } + k_kdot = 0.0; + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + //k_kdot += kdot[j][k] * Ku[j][k]; + k_kdot += Kdot[j][k] * Ku[j][k]; + } - beta_term = 0.0; - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) { - double D2beta = 0.0; - D2beta += d2beta[j][k][l]; + beta_term = 0.0; + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) { + double D2beta = 0.0; - for (int m = 0; m < 3; m++) { - D2beta += dG[j][l][k][m] * beta[m] + G[l][k][m] * dbeta[j][m] + G[l][j][m] * dbeta[k][m] - G[m][j][k] * dbeta[m][l]; + D2beta += d2beta[j][k][l]; - for (int n = 0; n < 3; n++) - D2beta += G[l][j][m] * G[m][k][n] * beta[n] - G[m][j][k] * G[l][m][n] * beta[n]; - } + for (int m = 0; m < 3; m++) { + D2beta += dG[j][l][k][m] * beta[m] + G[l][k][m] * dbeta[j][m] + G[l][j][m] * dbeta[k][m] - G[m][j][k] * dbeta[m][l]; - beta_term += -gu[j][k] * D2beta * dalpha[l]; + for (int n = 0; n < 3; n++) + D2beta += G[l][j][m] * G[m][k][n] * beta[n] - G[m][j][k] * G[l][m][n] * beta[n]; } - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - beta_term += -beta[k] * Ricm[j][k] * dalpha[j]; - - solver->diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_dc] = gu[0][0] + (on_axis ? gu[1][1] : 0.0); - solver->diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_dc] = gu[2][2]; - solver->diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_dc] = 2.0 * gu[0][2]; - solver->diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_dc] = -X[0] + (on_axis ? 0.0 : gu[1][1] / x); - solver->diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_dc] = -X[2]; - solver->diff_coeffs[MG2D_DIFF_COEFF_00]->data[idx_dc] = -k2; - - solver->rhs[idx_rhs] = -2.0 * alpha * Kij_Dij_alpha + Gammadot_term + - 2.0 * alpha * (k_kdot + 2.0 * alpha * k3) + beta_term; - - if (on_axis) { - solver->diff_coeffs[MG2D_DIFF_COEFF_20]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1]; - solver->diff_coeffs[MG2D_DIFF_COEFF_10]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1]; - } + + beta_term += -gu[j][k] * D2beta * dalpha[l]; + } + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + beta_term += -beta[k] * Ricm[j][k] * dalpha[j]; + + solver->diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_dc] = gu[0][0] + (on_axis ? gu[1][1] : 0.0); + solver->diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_dc] = gu[2][2]; + solver->diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_dc] = 2.0 * gu[0][2]; + solver->diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_dc] = -X[0] + (on_axis ? 0.0 : gu[1][1] / x); + solver->diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_dc] = -X[2]; + solver->diff_coeffs[MG2D_DIFF_COEFF_00]->data[idx_dc] = -k2; + + solver->rhs[idx_rhs] = -2.0 * alpha * Kij_Dij_alpha + Gammadot_term + + 2.0 * alpha * (k_kdot + 2.0 * alpha * k3) + beta_term; + + if (on_axis) { + solver->diff_coeffs[MG2D_DIFF_COEFF_20]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1]; + solver->diff_coeffs[MG2D_DIFF_COEFF_10]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1]; } + } } -- cgit v1.2.3