aboutsummaryrefslogtreecommitdiff
path: root/td_constraints.c
diff options
context:
space:
mode:
Diffstat (limited to 'td_constraints.c')
-rw-r--r--td_constraints.c195
1 files changed, 195 insertions, 0 deletions
diff --git a/td_constraints.c b/td_constraints.c
new file mode 100644
index 0000000..9ef6a80
--- /dev/null
+++ b/td_constraints.c
@@ -0,0 +1,195 @@
+/*
+ * Copyright 2017 Anton Khirnov <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <math.h>
+
+#include "common.h"
+#include "pssolve.h"
+#include "td_constraints.h"
+
+#define I_C sqrt(225.0 / (8.0 * M_PI))
+static double I_m2(double x, double a)
+{
+ return a * I_C * exp(-SQR(x));
+}
+static double I_m1(double x, double a)
+{
+ return a * I_C * (-2.0 * x * exp(-SQR(x)));
+}
+static double I_0(double x, double a)
+{
+ return a * I_C * (-2.0 + 4.0 * SQR(x)) * exp(-SQR(x));
+}
+static double I_p1(double x, double a)
+{
+ return a * I_C * (12.0 * x - 8.0 * x * SQR(x)) * exp(-SQR(x));
+}
+static double I_p2(double x, double a)
+{
+ return a * I_C * (12.0 - 48.0 * SQR(x) + 16.0 * SQR(SQR(x))) * exp(-SQR(x));
+}
+static double I_p3(double x, double a)
+{
+ return a * I_C * (-120.0 * x + 144.0 * x * SQR(x) - 32.0 * x * SQR(SQR(x))) * exp(-SQR(x));
+}
+
+static double K_rtheta(double a, double r, double theta)
+{
+ double r2 = r * r;
+ double r3 = r2 * r;
+
+ // equal to the below
+ return a * I_C * exp(-r2) * sin(2 * theta) * (32.0 * r3 - 48.0 * r);
+
+ //return ((I_p2(-r, a) + I_p2(r, a)) * r3 + 3.0 * (I_p1(-r, a) - I_p1(r, a)) * r2 +
+ // 6.0 * (I_0(-r, a) + I_0(r, a)) * r + 6.0 * (I_m1(-r, a) - I_m1(r, a))) * sin(2.0 * theta) / r4;
+}
+
+static double K_rtheta_dtheta(double a, double r, double theta)
+{
+ return 2.0 * cos(2.0 * theta) * K_rtheta(a, r, theta) / sin(2.0 * theta);
+}
+
+static double K_rtheta_dr(double a, double r, double theta)
+{
+ double r2 = r * r;
+ return -2 * r * K_rtheta(a, r, theta) + a * I_C * exp(-r2) * sin(2 * theta) * (3.0 * 32.0 * r2 - 48.0);
+ //return ((-I_p3(-r, a) + I_p3(r, a)) / r + 3.0 * (-I_p2(-r, a) - I_p2(r, a)) / (SQR(r)) +
+ // 6.0 * (-I_p1(-r, a) + I_p1(r, a)) / (r * SQR(r)) + 6.0 * (-I_0(-r, a) - I_0(r, a)) / (SQR(SQR(r))) +
+ // (I_p2(-r, a) + I_p2(r, a)) / (-SQR(r)) + 3.0 * (I_p1(-r, a) - I_p1(r, a)) * (-2.0 / (r * SQR(r))) +
+ // 6.0 * (I_0(-r, a) + I_0(r, a)) * (-3.0 / (SQR(SQR(r)))) + 6.0 * (I_m1(-r, a) - I_m1(r, a)) * (-4.0 / (r * SQR(SQR(r))))) * sin(2.0 * theta);
+}
+
+#include "eval_constraints.c"
+#include "eval_k_rtheta.c"
+
+//typedef void (*eval_metric_func)(double a, double L, double r_0, double x_0, double x_1, double psi_val_00, double psi_val_01, double psi_val_02, double psi_val_10, double psi_val_11, double psi_val_20, double *res);
+//
+//static const eval_metric_func eval_dmetric[3] = {
+// eval_dmetric_0,
+// eval_dmetric_1,
+// eval_dmetric_2,
+//};
+//
+//static const eval_metric_func eval_d2metric[3][3] = {
+// {
+// eval_d2metric_00,
+// eval_d2metric_01,
+// eval_d2metric_02,
+// },
+// {
+// eval_d2metric_10,
+// eval_d2metric_11,
+// eval_d2metric_12,
+// },
+// {
+// eval_d2metric_20,
+// eval_d2metric_21,
+// eval_d2metric_22,
+// },
+//};
+
+#define EQNAME constraint_eq_ham
+#define EQUATION TD_CONSTRAINT_EQ_HAM
+#include "td_constraints_template.c"
+#undef EQNAME
+#undef EQUATION
+
+#define EQNAME constraint_eq_mom_0
+#define EQUATION TD_CONSTRAINT_EQ_MOM_0
+#include "td_constraints_template.c"
+#undef EQNAME
+#undef EQUATION
+
+#define EQNAME constraint_eq_mom_1
+#define EQUATION TD_CONSTRAINT_EQ_MOM_1
+#include "td_constraints_template.c"
+#undef EQNAME
+#undef EQUATION
+
+void
+(*constraint_funcs[TD_CONSTRAINT_EQ_NB])(const double a,
+ unsigned int nb_coords[2], double *coords[2],
+ double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+ double *out) = {
+ [TD_CONSTRAINT_EQ_HAM] = constraint_eq_ham,
+ [TD_CONSTRAINT_EQ_MOM_0] = constraint_eq_mom_0,
+ [TD_CONSTRAINT_EQ_MOM_1] = constraint_eq_mom_1,
+};
+
+static double q_gundlach(const double a, double rho, double z)
+{
+ return a * SQR(rho) * exp(-(SQR(rho) + SQR(z)));
+}
+
+static double dq_rho_gundlach(const double a, double rho, double z)
+{
+ return a * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1.0 - rho * (rho));
+}
+
+static double d2q_rho_gundlach(const double a, double rho, double z)
+{
+ double rho2 = SQR(rho);
+ return a * 2 * exp(-SQR(rho) - SQR(z)) * (1.0 - 4.0 * rho * (rho) - rho2 + 2.0 * rho2 * SQR(rho));
+}
+
+static double d2q_rho_z_gundlach(const double a, double rho, double z)
+{
+ return a * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (rho * (rho) - 1.0);
+}
+
+static double dq_z_gundlach(const double a, double rho, double z)
+{
+ return -a * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z));
+}
+
+static double d2q_z_gundlach(const double a, double rho, double z)
+{
+ return a * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1.0);
+}
+
+
+int tdi_constraints_eval(enum TDConstraintEq eq,
+ const double a,
+ unsigned int nb_coords[2], double *coords[2],
+ double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+ double *out)
+{
+ //for (int idx1 = 0; idx1 < nb_coords[1]; idx1++) {
+ // const double theta = coords[1][idx1];
+
+ // for (int idx0 = 0; idx0 < nb_coords[0]; idx0++) {
+ // const double r = coords[0][idx0];
+ // const double rho = r * sin(theta);
+ // const double z = r * cos(theta);
+
+ // const double ct = cos(theta);
+ // const double st = sin(theta);
+
+ // const int idx = idx1 * nb_coords[0] + idx0;
+
+ // const double laplace = vars[eq][PSSOLVE_DIFF_ORDER_20][idx] +
+ // vars[eq][PSSOLVE_DIFF_ORDER_10][idx] * (2.0 / r) +
+ // vars[eq][PSSOLVE_DIFF_ORDER_02][idx] * (1.0 / SQR(r)) +
+ // vars[eq][PSSOLVE_DIFF_ORDER_01][idx] * (ct / (SQR(r) * st));
+
+ // out[idx] = laplace + 0.25 * (d2q_z_gundlach(a, rho, z) + d2q_rho_gundlach(a, rho, z)) * (1.0 + vars[eq][PSSOLVE_DIFF_ORDER_00][idx]);
+ // }
+ //}
+ constraint_funcs[eq](a, nb_coords, coords, vars, out);
+ return 0;
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-02 11:18:58 +0000