| 942 | } |
| 943 | |
| 944 | int Coordinates::geometry(bool recalculate_staggered, |
| 945 | bool force_interpolate_from_centre) { |
| 946 | TRACE("Coordinates::geometry"); |
| 947 | communicate(dx, dy, dz, g11, g22, g33, g12, g13, g23, g_11, g_22, g_33, g_12, g_13, |
| 948 | g_23, J, Bxy); |
| 949 | |
| 950 | output_progress.write("Calculating differential geometry terms\n"); |
| 951 | |
| 952 | if (min(abs(dx)) < 1e-8) { |
| 953 | throw BoutException("dx magnitude less than 1e-8"); |
| 954 | } |
| 955 | |
| 956 | if (min(abs(dy)) < 1e-8) { |
| 957 | throw BoutException("dy magnitude less than 1e-8"); |
| 958 | } |
| 959 | |
| 960 | if (min(abs(dz)) < 1e-8) { |
| 961 | throw BoutException("dz magnitude less than 1e-8"); |
| 962 | } |
| 963 | |
| 964 | // Check input metrics |
| 965 | checkContravariant(); |
| 966 | checkCovariant(); |
| 967 | |
| 968 | // Calculate Christoffel symbol terms (18 independent values) |
| 969 | // Note: This calculation is completely general: metric |
| 970 | // tensor can be 2D or 3D. For 2D, all DDZ terms are zero |
| 971 | |
| 972 | G1_11 = 0.5 * g11 * DDX(g_11) + g12 * (DDX(g_12) - 0.5 * DDY(g_11)) |
| 973 | + g13 * (DDX(g_13) - 0.5 * DDZ(g_11)); |
| 974 | G1_22 = g11 * (DDY(g_12) - 0.5 * DDX(g_22)) + 0.5 * g12 * DDY(g_22) |
| 975 | + g13 * (DDY(g_23) - 0.5 * DDZ(g_22)); |
| 976 | G1_33 = g11 * (DDZ(g_13) - 0.5 * DDX(g_33)) + g12 * (DDZ(g_23) - 0.5 * DDY(g_33)) |
| 977 | + 0.5 * g13 * DDZ(g_33); |
| 978 | G1_12 = 0.5 * g11 * DDY(g_11) + 0.5 * g12 * DDX(g_22) |
| 979 | + 0.5 * g13 * (DDY(g_13) + DDX(g_23) - DDZ(g_12)); |
| 980 | G1_13 = 0.5 * g11 * DDZ(g_11) + 0.5 * g12 * (DDZ(g_12) + DDX(g_23) - DDY(g_13)) |
| 981 | + 0.5 * g13 * DDX(g_33); |
| 982 | G1_23 = 0.5 * g11 * (DDZ(g_12) + DDY(g_13) - DDX(g_23)) |
| 983 | + 0.5 * g12 * (DDZ(g_22) + DDY(g_23) - DDY(g_23)) |
| 984 | // + 0.5 *g13*(DDZ(g_32) + DDY(g_33) - DDZ(g_23)); |
| 985 | // which equals |
| 986 | + 0.5 * g13 * DDY(g_33); |
| 987 | |
| 988 | G2_11 = 0.5 * g12 * DDX(g_11) + g22 * (DDX(g_12) - 0.5 * DDY(g_11)) |
| 989 | + g23 * (DDX(g_13) - 0.5 * DDZ(g_11)); |
| 990 | G2_22 = g12 * (DDY(g_12) - 0.5 * DDX(g_22)) + 0.5 * g22 * DDY(g_22) |
| 991 | + g23 * (DDY(g23) - 0.5 * DDZ(g_22)); |
| 992 | G2_33 = g12 * (DDZ(g_13) - 0.5 * DDX(g_33)) + g22 * (DDZ(g_23) - 0.5 * DDY(g_33)) |
| 993 | + 0.5 * g23 * DDZ(g_33); |
| 994 | G2_12 = 0.5 * g12 * DDY(g_11) + 0.5 * g22 * DDX(g_22) |
| 995 | + 0.5 * g23 * (DDY(g_13) + DDX(g_23) - DDZ(g_12)); |
| 996 | G2_13 = |
| 997 | // 0.5 *g21*(DDZ(g_11) + DDX(g_13) - DDX(g_13)) |
| 998 | // which equals |
| 999 | 0.5 * g12 * (DDZ(g_11) + DDX(g_13) - DDX(g_13)) |
| 1000 | // + 0.5 *g22*(DDZ(g_21) + DDX(g_23) - DDY(g_13)) |
| 1001 | // which equals |