From 8a3cb970c5c437b471ae32f8c45762ec3012ae50 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 5 Oct 2018 15:38:21 +0200 Subject: td_constraints: tie the constraint equations to the family Preparation for the ingoing pulse family, which has a slightly different metric expression and so different form of the constraint equations. --- td_constraints.c | 159 +++++++++++++++++++++++++++++-------------------------- 1 file changed, 83 insertions(+), 76 deletions(-) diff --git a/td_constraints.c b/td_constraints.c index 41c1c5b..5618cb3 100644 --- a/td_constraints.c +++ b/td_constraints.c @@ -27,79 +27,19 @@ struct TDConstraintEvalPriv { double *dk_rtheta[2]; }; +typedef void (*ConstraintEvalCB)(TDConstraintEvalContext *ctx, + double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], + double *out); + typedef struct TDFamilyDef { double (*eval_krt) (TDConstraintEvalContext *ctx, double r, double theta); double (*eval_krt_dr)(TDConstraintEvalContext *ctx, double r, double theta); double (*eval_krt_dt)(TDConstraintEvalContext *ctx, double r, double theta); double a_converge; double a_diverge; -} TDFamilyDef; - -static double k_rtheta_eval_ae_time_antisym(TDConstraintEvalContext *ctx, - double r, double theta) -{ - const double r2 = SQR(r); - const double r3 = r * r2; - return -60.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI); -} - -static double k_rtheta_eval_dr_ae_time_antisym(TDConstraintEvalContext *ctx, - double r, double theta) -{ - const double r2 = SQR(r); - const double r4 = SQR(r2); - return -60.0 * sqrt(2) * (3.0 - 12.0 * r2 + 4.0 * r4) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI); -} - -static double k_rtheta_eval_dt_ae_time_antisym(TDConstraintEvalContext *ctx, - double r, double theta) -{ - const double r2 = SQR(r); - const double r3 = r * r2; - return -120.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * cos(2.0 * theta) * ctx->amplitude / sqrt(M_PI); -} - -static const TDFamilyDef ae_time_antisym = { - .eval_krt = k_rtheta_eval_ae_time_antisym, - .eval_krt_dr = k_rtheta_eval_dr_ae_time_antisym, - .eval_krt_dt = k_rtheta_eval_dt_ae_time_antisym, - .a_converge = 0.0092078837420, - .a_diverge = 0.0092078837427, -}; - -static double k_rtheta_eval_simple_time_antisym(TDConstraintEvalContext *ctx, - double r, double theta) -{ - const double r2 = SQR(r); - return ctx->amplitude * r * exp(-r2) * sin(2.0 * theta); -} - -static double k_rtheta_eval_dr_simple_time_antisym(TDConstraintEvalContext *ctx, - double r, double theta) -{ - const double r2 = SQR(r); - return ctx->amplitude * (1.0 - 2.0 * r2) * exp(-r2) * sin(2.0 * theta); -} - -static double k_rtheta_eval_dt_simple_time_antisym(TDConstraintEvalContext *ctx, - double r, double theta) -{ - const double r2 = SQR(r); - return 2.0 * ctx->amplitude * r * exp(-r2) * cos(2.0 * theta); -} -static const TDFamilyDef simple_time_antisym = { - .eval_krt = k_rtheta_eval_simple_time_antisym, - .eval_krt_dr = k_rtheta_eval_dr_simple_time_antisym, - .eval_krt_dt = k_rtheta_eval_dt_simple_time_antisym, - .a_converge = 1.018918628476058, - .a_diverge = 1.018918628476968, -}; - -static const TDFamilyDef *td_families[] = { - [TD_FAMILY_AE_TIME_ANTISYM] = &ae_time_antisym, - [TD_FAMILY_SIMPLE_TIME_ANTISYM] = &simple_time_antisym, -}; + const ConstraintEvalCB *constraint_eval; +} TDFamilyDef; static void constraint_eval_ham_confflat(TDConstraintEvalContext *ctx, @@ -212,6 +152,81 @@ constraint_eval_mom_t_confflat(TDConstraintEvalContext *ctx, } } +static const ConstraintEvalCB constraint_funcs_confflat[TD_CONSTRAINT_EQ_NB] = { + [TD_CONSTRAINT_EQ_HAM] = constraint_eval_ham_confflat, + [TD_CONSTRAINT_EQ_MOM_0] = constraint_eval_mom_r_confflat, + [TD_CONSTRAINT_EQ_MOM_1] = constraint_eval_mom_t_confflat, +}; + +static double k_rtheta_eval_ae_time_antisym(TDConstraintEvalContext *ctx, + double r, double theta) +{ + const double r2 = SQR(r); + const double r3 = r * r2; + return -60.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI); +} + +static double k_rtheta_eval_dr_ae_time_antisym(TDConstraintEvalContext *ctx, + double r, double theta) +{ + const double r2 = SQR(r); + const double r4 = SQR(r2); + return -60.0 * sqrt(2) * (3.0 - 12.0 * r2 + 4.0 * r4) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI); +} + +static double k_rtheta_eval_dt_ae_time_antisym(TDConstraintEvalContext *ctx, + double r, double theta) +{ + const double r2 = SQR(r); + const double r3 = r * r2; + return -120.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * cos(2.0 * theta) * ctx->amplitude / sqrt(M_PI); +} + +static const TDFamilyDef ae_time_antisym = { + .eval_krt = k_rtheta_eval_ae_time_antisym, + .eval_krt_dr = k_rtheta_eval_dr_ae_time_antisym, + .eval_krt_dt = k_rtheta_eval_dt_ae_time_antisym, + .a_converge = 0.0092078837420, + .a_diverge = 0.0092078837427, + .constraint_eval = constraint_funcs_confflat, +}; + +static double k_rtheta_eval_simple_time_antisym(TDConstraintEvalContext *ctx, + double r, double theta) +{ + const double r2 = SQR(r); + return ctx->amplitude * r * exp(-r2) * sin(2.0 * theta); +} + +static double k_rtheta_eval_dr_simple_time_antisym(TDConstraintEvalContext *ctx, + double r, double theta) +{ + const double r2 = SQR(r); + return ctx->amplitude * (1.0 - 2.0 * r2) * exp(-r2) * sin(2.0 * theta); +} + +static double k_rtheta_eval_dt_simple_time_antisym(TDConstraintEvalContext *ctx, + double r, double theta) +{ + const double r2 = SQR(r); + return 2.0 * ctx->amplitude * r * exp(-r2) * cos(2.0 * theta); +} + +static const TDFamilyDef simple_time_antisym = { + .eval_krt = k_rtheta_eval_simple_time_antisym, + .eval_krt_dr = k_rtheta_eval_dr_simple_time_antisym, + .eval_krt_dt = k_rtheta_eval_dt_simple_time_antisym, + //.a_converge = 1.018918628476058, + .a_converge = 1.0189186, + .a_diverge = 1.018918628476968, + .constraint_eval = constraint_funcs_confflat, +}; + +static const TDFamilyDef *td_families[] = { + [TD_FAMILY_AE_TIME_ANTISYM] = &ae_time_antisym, + [TD_FAMILY_SIMPLE_TIME_ANTISYM] = &simple_time_antisym, +}; + double tdi_constraint_eval_k_rtheta(TDConstraintEvalContext *ctx, double r, double theta) { const TDFamilyDef *fd = td_families[ctx->family]; @@ -310,19 +325,11 @@ void tdi_constraint_eval_free(TDConstraintEvalContext **pctx) *pctx = NULL; } -static const void -(*constraint_funcs[TD_CONSTRAINT_EQ_NB])(TDConstraintEvalContext *ctx, - double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], - double *out) = { - [TD_CONSTRAINT_EQ_HAM] = constraint_eval_ham_confflat, - [TD_CONSTRAINT_EQ_MOM_0] = constraint_eval_mom_r_confflat, - [TD_CONSTRAINT_EQ_MOM_1] = constraint_eval_mom_t_confflat, -}; - int tdi_constraint_eval(TDConstraintEvalContext *ctx, enum TDConstraintEq eq, const double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], double *out) { - constraint_funcs[eq](ctx, vars, out); + const TDFamilyDef *fd = td_families[ctx->family]; + fd->constraint_eval[eq](ctx, vars, out); return 0; } -- cgit v1.2.3