From 222dc8ae57c7eba0967f2f932b78961cafb0ee83 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Wed, 23 Nov 2016 17:21:38 +0100 Subject: solve: add the shift terms --- src/qms_solve.c | 591 ++++++++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 537 insertions(+), 54 deletions(-) diff --git a/src/qms_solve.c b/src/qms_solve.c index 7136f38..c84c1d8 100644 --- a/src/qms_solve.c +++ b/src/qms_solve.c @@ -73,7 +73,6 @@ enum MetricVars { XTDOT_X, XTDOT_Y, XTDOT_Z, - PHIDOT, NB_METRIC_VARS, }; @@ -85,23 +84,66 @@ enum InterpMetricVars { I_GTXY, I_GTXZ, I_GTYZ, + I_GTXX_DX, + I_GTYY_DX, + I_GTZZ_DX, + I_GTXZ_DX, + I_GTXX_DZ, + I_GTYY_DZ, + I_GTZZ_DZ, + I_GTXZ_DZ, + I_GTXX_DXX, + I_GTYY_DXX, + I_GTZZ_DXX, + I_GTXZ_DXX, + I_GTXX_DXZ, + I_GTYY_DXZ, + I_GTZZ_DXZ, + I_GTXZ_DXZ, + I_GTXX_DZZ, + I_GTYY_DZZ, + I_GTZZ_DZZ, + I_GTXZ_DZZ, I_PHI, I_PHI_DX, I_PHI_DY, I_PHI_DZ, + I_PHI_DXX, + I_PHI_DZZ, + I_PHI_DXZ, I_ATXX, I_ATYY, I_ATZZ, I_ATXY, I_ATXZ, I_ATYZ, + I_ATXX_DX, + I_ATYY_DX, + I_ATZZ_DX, + I_ATXZ_DX, + I_ATXX_DZ, + I_ATYY_DZ, + I_ATZZ_DZ, + I_ATXZ_DZ, I_K, + I_TRK_DX, + I_TRK_DZ, I_XTX, I_XTY, I_XTZ, I_BETAX, I_BETAY, I_BETAZ, + I_BETAX_DX, + I_BETAX_DZ, + I_BETAX_DXX, + I_BETAX_DXZ, + I_BETAX_DZZ, + I_BETAZ_DX, + I_BETAZ_DZ, + I_BETAZ_DXX, + I_BETAZ_DXZ, + I_BETAZ_DZZ, I_ALPHA, I_ALPHA_DX, I_ALPHA_DY, @@ -121,10 +163,6 @@ enum InterpMetricVars { I_XTDOT_X, I_XTDOT_Y, I_XTDOT_Z, - I_PHIDOT, - I_PHIDOT_DX, - I_PHIDOT_DY, - I_PHIDOT_DZ, NB_INTERP_VARS, }; @@ -199,7 +237,6 @@ static const char *metric_vars[] = { [XTDOT_X] = "ML_CCZ4::Xtdot1", [XTDOT_Y] = "ML_CCZ4::Xtdot2", [XTDOT_Z] = "ML_CCZ4::Xtdot3", - [PHIDOT] = "ML_CCZ4::phidot", #else [GTXX] = "ML_BSSN::gt11", [GTYY] = "ML_BSSN::gt22", @@ -232,7 +269,6 @@ static const char *metric_vars[] = { [XTDOT_X] = "ML_BSSN::Xtdot1", [XTDOT_Y] = "ML_BSSN::Xtdot2", [XTDOT_Z] = "ML_BSSN::Xtdot3", - [PHIDOT] = "ML_BSSN::phidot", #endif }; @@ -244,23 +280,66 @@ static const CCTK_INT interp_operation_indices[] = { [I_GTXY] = GTXY, [I_GTXZ] = GTXZ, [I_GTYZ] = GTYZ, + [I_GTXX_DX] = GTXX, + [I_GTYY_DX] = GTYY, + [I_GTZZ_DX] = GTZZ, + [I_GTXZ_DX] = GTXZ, + [I_GTXX_DZ] = GTXX, + [I_GTYY_DZ] = GTYY, + [I_GTZZ_DZ] = GTZZ, + [I_GTXZ_DZ] = GTXZ, + [I_GTXX_DXX] = GTXX, + [I_GTYY_DXX] = GTYY, + [I_GTZZ_DXX] = GTZZ, + [I_GTXZ_DXX] = GTXZ, + [I_GTXX_DXZ] = GTXX, + [I_GTYY_DXZ] = GTYY, + [I_GTZZ_DXZ] = GTZZ, + [I_GTXZ_DXZ] = GTXZ, + [I_GTXX_DZZ] = GTXX, + [I_GTYY_DZZ] = GTYY, + [I_GTZZ_DZZ] = GTZZ, + [I_GTXZ_DZZ] = GTXZ, [I_PHI] = PHI, [I_PHI_DX] = PHI, [I_PHI_DY] = PHI, [I_PHI_DZ] = PHI, + [I_PHI_DXX] = PHI, + [I_PHI_DZZ] = PHI, + [I_PHI_DXZ] = PHI, [I_ATXX] = ATXX, [I_ATYY] = ATYY, [I_ATZZ] = ATZZ, [I_ATXY] = ATXY, [I_ATXZ] = ATXZ, [I_ATYZ] = ATYZ, + [I_ATXX_DX] = ATXX, + [I_ATYY_DX] = ATYY, + [I_ATZZ_DX] = ATZZ, + [I_ATXZ_DX] = ATXZ, + [I_ATXX_DZ] = ATXX, + [I_ATYY_DZ] = ATYY, + [I_ATZZ_DZ] = ATZZ, + [I_ATXZ_DZ] = ATXZ, [I_K] = K, + [I_TRK_DX] = K, + [I_TRK_DZ] = K, [I_XTX] = XTX, [I_XTY] = XTY, [I_XTZ] = XTZ, - [I_BETAX] = BETAX, - [I_BETAY] = BETAY, - [I_BETAZ] = BETAZ, + [I_BETAX] = BETAX, + [I_BETAY] = BETAY, + [I_BETAZ] = BETAZ, + [I_BETAX_DX] = BETAX, + [I_BETAX_DZ] = BETAX, + [I_BETAX_DXX] = BETAX, + [I_BETAX_DXZ] = BETAX, + [I_BETAX_DZZ] = BETAX, + [I_BETAZ_DX] = BETAZ, + [I_BETAZ_DZ] = BETAZ, + [I_BETAZ_DXX] = BETAZ, + [I_BETAZ_DXZ] = BETAZ, + [I_BETAZ_DZZ] = BETAZ, [I_ALPHA] = ALPHA, [I_ALPHA_DX] = ALPHA, [I_ALPHA_DY] = ALPHA, @@ -280,10 +359,6 @@ static const CCTK_INT interp_operation_indices[] = { [I_XTDOT_X] = XTDOT_X, [I_XTDOT_Y] = XTDOT_Y, [I_XTDOT_Z] = XTDOT_Z, - [I_PHIDOT] = PHIDOT, - [I_PHIDOT_DX] = PHIDOT, - [I_PHIDOT_DY] = PHIDOT, - [I_PHIDOT_DZ] = PHIDOT, }; /* the operation (plain value or x/y/z-derivative) to apply during interpolation */ @@ -294,23 +369,66 @@ static const CCTK_INT interp_operation_codes[] = { [I_GTXY] = 0, [I_GTXZ] = 0, [I_GTYZ] = 0, + [I_GTXX_DX] = 1, + [I_GTYY_DX] = 1, + [I_GTZZ_DX] = 1, + [I_GTXZ_DX] = 1, + [I_GTXX_DZ] = 3, + [I_GTYY_DZ] = 3, + [I_GTZZ_DZ] = 3, + [I_GTXZ_DZ] = 3, + [I_GTXX_DXX] = 11, + [I_GTYY_DXX] = 11, + [I_GTZZ_DXX] = 11, + [I_GTXZ_DXX] = 11, + [I_GTXX_DXZ] = 13, + [I_GTYY_DXZ] = 13, + [I_GTZZ_DXZ] = 13, + [I_GTXZ_DXZ] = 13, + [I_GTXX_DZZ] = 33, + [I_GTYY_DZZ] = 33, + [I_GTZZ_DZZ] = 33, + [I_GTXZ_DZZ] = 33, [I_PHI] = 0, [I_PHI_DX] = 1, [I_PHI_DY] = 2, [I_PHI_DZ] = 3, + [I_PHI_DXX] = 11, + [I_PHI_DZZ] = 33, + [I_PHI_DXZ] = 13, [I_ATXX] = 0, [I_ATYY] = 0, [I_ATZZ] = 0, [I_ATXY] = 0, [I_ATXZ] = 0, [I_ATYZ] = 0, + [I_ATXX_DX] = 1, + [I_ATYY_DX] = 1, + [I_ATZZ_DX] = 1, + [I_ATXZ_DX] = 1, + [I_ATXX_DZ] = 3, + [I_ATYY_DZ] = 3, + [I_ATZZ_DZ] = 3, + [I_ATXZ_DZ] = 3, [I_K] = 0, + [I_TRK_DX] = 1, + [I_TRK_DZ] = 3, [I_XTX] = 0, [I_XTY] = 0, [I_XTZ] = 0, - [I_BETAX] = 0, - [I_BETAY] = 0, - [I_BETAZ] = 0, + [I_BETAX] = 0, + [I_BETAY] = 0, + [I_BETAZ] = 0, + [I_BETAX_DX] = 1, + [I_BETAX_DZ] = 3, + [I_BETAX_DXX] = 11, + [I_BETAX_DXZ] = 13, + [I_BETAX_DZZ] = 33, + [I_BETAZ_DX] = 1, + [I_BETAZ_DZ] = 3, + [I_BETAZ_DXX] = 11, + [I_BETAZ_DXZ] = 13, + [I_BETAZ_DZZ] = 33, [I_ALPHA] = 0, [I_ALPHA_DX] = 1, [I_ALPHA_DY] = 2, @@ -330,10 +448,6 @@ static const CCTK_INT interp_operation_codes[] = { [I_XTDOT_X] = 0, [I_XTDOT_Y] = 0, [I_XTDOT_Z] = 0, - [I_PHIDOT] = 0, - [I_PHIDOT_DX] = 1, - [I_PHIDOT_DY] = 2, - [I_PHIDOT_DZ] = 3, }; /* interpolate the cactus gridfunctions onto the pseudospectral grid */ @@ -363,8 +477,16 @@ static int calc_eq_coeffs(QMSSolver *ctx) const double z = s->interp_coords[2][i]; const int zaxis = x <= EPS; - double Am[3][3], K[3][3], Km[3][3], Ku[3][3], gtu[3][3]; - double k2, kij_dij_alpha, k_kdot, k3; + 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]; + double Ric[3][3], Ricm[3][3]; + double dgdot[3][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]; const double gtxx = s->interp_values[I_GTXX][i]; const double gtyy = s->interp_values[I_GTYY][i]; @@ -377,6 +499,111 @@ static int calc_eq_coeffs(QMSSolver *ctx) { gtxy, gtyy, gtyz }, { gtxz, gtyz, gtzz }}; + const double dx_gt11 = s->interp_values[I_GTXX_DX][i]; + const double dx_gt22 = s->interp_values[I_GTYY_DX][i]; + const double dx_gt33 = s->interp_values[I_GTZZ_DX][i]; + const double dx_gt13 = s->interp_values[I_GTXZ_DX][i]; + + const double dz_gt11 = s->interp_values[I_GTXX_DZ][i]; + const double dz_gt22 = s->interp_values[I_GTYY_DZ][i]; + const double dz_gt33 = s->interp_values[I_GTZZ_DZ][i]; + const double dz_gt13 = s->interp_values[I_GTXZ_DZ][i]; + + 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, zaxis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0 }, + { zaxis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0, zaxis ? dx_gt13 : gtxz / x }, + { 0.0, zaxis ? 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 = s->interp_values[I_GTXX_DXX][i]; + const double dxx_gt22 = s->interp_values[I_GTYY_DXX][i]; + const double dxx_gt33 = s->interp_values[I_GTZZ_DXX][i]; + const double dxx_gt13 = s->interp_values[I_GTXZ_DXX][i]; + + const double dxz_gt11 = s->interp_values[I_GTXX_DXZ][i]; + const double dxz_gt22 = s->interp_values[I_GTYY_DXZ][i]; + const double dxz_gt33 = s->interp_values[I_GTZZ_DXZ][i]; + const double dxz_gt13 = s->interp_values[I_GTXZ_DXZ][i]; + + const double dzz_gt11 = s->interp_values[I_GTXX_DZZ][i]; + const double dzz_gt22 = s->interp_values[I_GTYY_DZZ][i]; + const double dzz_gt33 = s->interp_values[I_GTZZ_DZZ][i]; + const double dzz_gt13 = s->interp_values[I_GTXZ_DZZ][i]; + + 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, zaxis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 }, + { zaxis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0, + zaxis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, + { 0.0, zaxis ? 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, zaxis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 }, + { zaxis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0, + zaxis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, + { 0.0, zaxis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 }, + }, + { + { zaxis ? dxx_gt22 : dx_gt11 / x - 2 * (gtxx - gtyy) / SQR(x), 0.0, + zaxis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, + { 0.0, zaxis ? dxx_gt11 : dx_gt22 / x + 2.0 * (gtxx - gtyy) / SQR(x), 0.0 }, + { zaxis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0, zaxis ? dxx_gt33 : dx_gt33 / x }, + }, + { + { 0.0, zaxis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 }, + { zaxis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0, + zaxis ? dxz_gt13 : dz_gt13 / x }, + { 0.0, zaxis ? dxz_gt13 : dz_gt13 / x, 0.0 }, + }, + + }, + { + { + { dxz_gt11, 0.0, dxz_gt13 }, + { 0.0, dxz_gt22, 0.0 }, + { dxz_gt13, 0.0, dxz_gt33 }, + }, + { + { 0.0, zaxis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 }, + { zaxis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0, + zaxis ? dxz_gt13 : dz_gt13 / x }, + { 0.0, zaxis ? 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 }, + }, + + }, + }; + const double Atxx = s->interp_values[I_ATXX][i]; const double Atyy = s->interp_values[I_ATYY][i]; const double Atzz = s->interp_values[I_ATZZ][i]; @@ -384,17 +611,62 @@ static int calc_eq_coeffs(QMSSolver *ctx) const double Atxz = s->interp_values[I_ATXZ][i]; const double Atyz = s->interp_values[I_ATYZ][i]; - const double phi = s->interp_values[I_PHI][i]; - - const double phidot = s->interp_values[I_PHIDOT][i]; - const double phidot_dx = s->interp_values[I_PHIDOT_DX][i]; - const double phidot_dz = s->interp_values[I_PHIDOT_DZ][i]; - const double At[3][3] = {{ Atxx, Atxy, Atxz }, { Atxy, Atyy, Atyz }, { Atxz, Atyz, Atzz }}; + const double dx_At11 = s->interp_values[I_ATXX_DX][i]; + const double dx_At22 = s->interp_values[I_ATYY_DX][i]; + const double dx_At33 = s->interp_values[I_ATZZ_DX][i]; + const double dx_At13 = s->interp_values[I_ATXZ_DX][i]; + + const double dz_At11 = s->interp_values[I_ATXX_DZ][i]; + const double dz_At22 = s->interp_values[I_ATYY_DZ][i]; + const double dz_At33 = s->interp_values[I_ATZZ_DZ][i]; + const double dz_At13 = s->interp_values[I_ATXZ_DZ][i]; + + 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, zaxis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0 }, + { zaxis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0, zaxis ? dx_At13 : Atxz / x }, + { 0.0, zaxis ? 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 = s->interp_values[I_PHI][i]; + + const double dx_phi = s->interp_values[I_PHI_DX][i]; + const double dz_phi = s->interp_values[I_PHI_DZ][i]; + + const double dphi[3] = { dx_phi, 0.0, dz_phi }; + + const double phi_dxx = s->interp_values[I_PHI_DXX][i]; + const double phi_dzz = s->interp_values[I_PHI_DZZ][i]; + const double phi_dxz = s->interp_values[I_PHI_DXZ][i]; + + const double d2phi[3][3] = { + { phi_dxx, 0.0, phi_dxz }, + { 0.0, zaxis ? phi_dxx : dx_phi / x, 0.0 }, + { phi_dxz, 0.0, phi_dzz }, + }; + const double trK = s->interp_values[I_K][i]; + + const double dx_trK = s->interp_values[I_TRK_DX][i]; + const double dz_trK = s->interp_values[I_TRK_DZ][i]; + + const double dtrK[3] = { dx_trK, 0.0, dz_trK }; + const double kdot_xx = s->interp_values[I_KDOT_XX][i]; const double kdot_yy = s->interp_values[I_KDOT_YY][i]; const double kdot_zz = s->interp_values[I_KDOT_ZZ][i]; @@ -409,13 +681,57 @@ static int calc_eq_coeffs(QMSSolver *ctx) const double alpha = s->interp_values[I_ALPHA][i]; const double dx_alpha = s->interp_values[I_ALPHA_DX][i]; const double dz_alpha = s->interp_values[I_ALPHA_DZ][i]; + + const double dalpha[3] = { dx_alpha, 0.0, dz_alpha }; + const double dxx_alpha = s->interp_values[I_ALPHA_DXX][i]; const double dzz_alpha = s->interp_values[I_ALPHA_DZZ][i]; const double dxz_alpha = s->interp_values[I_ALPHA_DXZ][i]; - const double dij_alpha[3][3] = {{ dxx_alpha, 0, dxz_alpha }, - { 0, zaxis ? dxx_alpha : dx_alpha / x, 0 }, - { dxz_alpha, 0, dzz_alpha }}; + const double d2alpha[3][3] = {{ dxx_alpha, 0, dxz_alpha }, + { 0, zaxis ? dxx_alpha : dx_alpha / x, 0 }, + { dxz_alpha, 0, dzz_alpha }}; + + const double betax = s->interp_values[I_BETAX][i]; + const double betaz = s->interp_values[I_BETAZ][i]; + + const double beta[3] = { betax, 0.0, betaz }; + + const double dx_betax = s->interp_values[I_BETAX_DX][i]; + const double dz_betax = s->interp_values[I_BETAX_DZ][i]; + + const double dx_betaz = s->interp_values[I_BETAZ_DX][i]; + const double dz_betaz = s->interp_values[I_BETAZ_DZ][i]; + + const double dbeta[3][3] = {{ dx_betax, 0.0, dx_betaz }, + { 0.0, zaxis ? dx_betax : betax / x, 0.0 }, + { dz_betax, 0.0, dz_betaz }}; + + const double dxx_betax = s->interp_values[I_BETAX_DXX][i]; + const double dxz_betax = s->interp_values[I_BETAX_DXZ][i]; + const double dzz_betax = s->interp_values[I_BETAX_DZZ][i]; + + const double dxx_betaz = s->interp_values[I_BETAZ_DXX][i]; + const double dxz_betaz = s->interp_values[I_BETAZ_DXZ][i]; + const double dzz_betaz = s->interp_values[I_BETAZ_DZZ][i]; + + const double d2beta[3][3][3] = { + { + { dxx_betax, 0.0, dxx_betaz }, + { 0.0, zaxis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 }, + { dxz_betax, 0.0, dxz_betaz }, + }, + { + { 0.0, zaxis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 }, + { zaxis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0, zaxis ? dxx_betaz : dx_betaz / x }, + { 0.0, zaxis ? dxz_betax : dz_betax / x, 0.0 }, + }, + { + { dxz_betax, 0.0, dxz_betaz }, + { 0.0, zaxis ? dxz_betax : dz_betax / x, 0.0 }, + { dzz_betax, 0.0, dzz_betaz }, + }, + }; const double Xtx = s->interp_values[I_XTX][i]; const double Xtz = s->interp_values[I_XTZ][i]; @@ -436,16 +752,107 @@ static int calc_eq_coeffs(QMSSolver *ctx) 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; + } + + // ∂_j γ_{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 γ^{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} + 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; + } + + // ∂_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 β_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; + } + + // ∂_{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 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; + } + + // K_{ij} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) - K[j][k] = At[j][k] / SQR(phi) + gt[j][k] * trK; + K[j][k] = At[j][k] / SQR(phi) + (1.0 / 3.0) * g[j][k] * trK; 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 += SQR(phi) * gtu[j][l] * K[l][k]; + val += gu[j][l] * K[l][k]; Km[j][k] = val; } @@ -454,23 +861,65 @@ static int calc_eq_coeffs(QMSSolver *ctx) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) - val += SQR(phi) * gtu[j][l] * Km[k][l]; + val += gu[j][l] * Km[k][l]; Ku[j][k] = val; } - // \tilde{A}_{i}^j + // ∂_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++) { - double val = 0.0; - for (int l = 0; l < 3; l++) - val += gtu[j][l] * At[l][k]; - Am[j][k] = val; + gudot[j][k] = 2.0 * alpha * Ku[j][k]; + + for (int l = 0; l < 3; l++) { + gudot[j][k] += -gu[j][l] * dbeta[l][k] - gu[k][l] * dbeta[l][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]; + } } - kij_dij_alpha = 0.0; + // \dot{Γ}^j_{kl} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) - kij_dij_alpha += Ku[j][k] * dij_alpha[j][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]); + + 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]; + } + + Kij_Dij_alpha = 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 += G[l][j][k] * dalpha[l]; + + Kij_Dij_alpha += Ku[j][k] * (d2alpha[j][k] - val); + } k_kdot = 0.0; for (int j = 0; j < 3; j++) @@ -482,8 +931,8 @@ static int calc_eq_coeffs(QMSSolver *ctx) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) - val += Km[k][l] * Ku[l][j]; - k3 += val * K[j][k]; + val += Km[k][l] * Km[l][j]; + k3 += val * Km[j][k]; } // K_{ij} K^{ij} @@ -492,6 +941,48 @@ static int calc_eq_coeffs(QMSSolver *ctx) 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 += G[l][l][m] * G[m][j][k] - G[l][k][m] * G[m][j][l]; + Ric[j][k] = 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; + } + + 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]; + + 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]; + + 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]; + } + + 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]; + { const double gtuxx = gtu[0][0]; const double gtuyy = gtu[1][1]; @@ -504,17 +995,9 @@ static int calc_eq_coeffs(QMSSolver *ctx) const double Xtx = s->interp_values[I_XTX][i]; const double Xtz = s->interp_values[I_XTZ][i]; - const double betax = s->interp_values[I_BETAX][i]; - const double betaz = s->interp_values[I_BETAZ][i]; - const double Xx = SQR(phi) * (Xtx + (phi_dx * gtuxx + phi_dz * gtuxz) / phi); const double Xz = SQR(phi) * (Xtz + (phi_dx * gtuxz + phi_dz * gtuzz) / phi); - const double Xdot_x = 2 * phi * phidot * Xtx + SQR(phi) * Xtdot_x + phi * (phidot_dx * gtuxx + phidot_dz * gtuxz) - - phidot * (phi_dx * gtuxx + phi_dz * gtuxz) + 2 * alpha * (phi_dx * Ku[0][0] + phi_dz * Ku[0][2]) / phi; - const double Xdot_z = 2 * phi * phidot * Xtz + SQR(phi) * Xtdot_z + phi * (phidot_dz * gtuzz + phidot_dx * gtuxz) - - phidot * (phi_dz * gtuzz + phi_dx * gtuxz) + 2 * alpha * (phi_dz * Ku[2][2] + phi_dx * Ku[0][2]) / phi; - s->eq_coeffs[PSSOLVE_DIFF_ORDER_20][i] = SQR(phi) * (gtuxx + ((x <= EPS) ? gtuyy : 0.0)); s->eq_coeffs[PSSOLVE_DIFF_ORDER_02][i] = SQR(phi) * gtuzz; s->eq_coeffs[PSSOLVE_DIFF_ORDER_11][i] = SQR(phi) * gtuxz * 2; @@ -522,8 +1005,8 @@ static int calc_eq_coeffs(QMSSolver *ctx) s->eq_coeffs[PSSOLVE_DIFF_ORDER_01][i] = -Xz; s->eq_coeffs[PSSOLVE_DIFF_ORDER_00][i] = -k2; - s->rhs[i] = -2 * alpha * kij_dij_alpha + Xdot_x * dx_alpha + Xdot_z * dz_alpha + - 2 * (k_kdot + 2 * alpha * k3) * alpha; + s->rhs[i] = -2 * alpha * Kij_Dij_alpha + Gammadot_term + + 2 * alpha * (k_kdot + 2 * alpha * k3) + beta_term; } } -- cgit v1.2.3