summaryrefslogtreecommitdiff
path: root/src/qms.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/qms.c')
-rw-r--r--src/qms.c52
1 files changed, 37 insertions, 15 deletions
diff --git a/src/qms.c b/src/qms.c
index 95c645c..b5b8100 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -21,6 +21,43 @@
#include "qms.h"
#include "qms_solve.h"
+/* precomputed values for a given refined grid */
+typedef struct CoordPatch {
+ CCTK_REAL origin[3];
+ CCTK_INT delta[3];
+ CCTK_INT size[3];
+
+ // basis values on the grid
+ double *basis_val_r;
+ double *basis_val_z;
+
+ double *transform_z;
+ double *transform_matrix;
+ double *transform_matrix1;
+ double *transform_tmp;
+
+ int y_idx;
+} CoordPatch;
+
+struct QMSContext {
+ QMSSolver *solver;
+ cGH *gh;
+
+ struct {
+ double time;
+ double *coeffs;
+ } solution_cache[8];
+ int nb_solutions;
+
+ double *coeffs_eval;
+
+ uint64_t grid_expand_count;
+ uint64_t grid_expand_time;
+
+ CoordPatch *patches;
+ int nb_patches;
+};
+
double scale_factor;
/* get an approximate "main" frequency component in a basis function */
@@ -37,7 +74,6 @@ static CoordPatch *get_coord_patch(QMSContext *ms,
CoordPatch *cp;
int64_t grid_size;
- int i;
for (int i = 0; i < ms->nb_patches; i++) {
cp = &ms->patches[i];
@@ -138,20 +174,6 @@ static CoordPatch *get_coord_patch(QMSContext *ms,
//cp->basis_val_z [i + ms->gh->cctk_lsh[2] * j] = ms->basis->eval(zz, j);
}
posix_memalign((void**)&cp->transform_z, 32, sizeof(*cp->transform_z) * cctkGH->cctk_lsh[2] * ms->solver->nb_coeffs[0]);
- posix_memalign((void**)&cp->one, 32, sizeof(*cp->one) * grid_size);
- for (int i = 0; i < grid_size; i++)
- cp->one[i] = 1.0;
-
- posix_memalign((void**)&cp->w_scale, 32, sizeof(*cp->w_scale) * grid_size);
- for (int k = 0; k < ms->gh->cctk_lsh[2]; k++)
- for (int j = 0; j < ms->gh->cctk_lsh[1]; j++)
- for (int i = 0; i < ms->gh->cctk_lsh[0]; i++) {
- int idx = CCTK_GFINDEX3D(ms->gh, i, j, k);
- double r = sqrt(SQR(x[idx]) + SQR(y[idx]) + SQR(z[idx]));
- const double R = 32.0;
- const double width = 4.0;
- cp->w_scale[idx] = (r > R) ? exp(-pow((r - R) / width, 4.0)) : 1.0;
- }
#endif
ms->nb_patches++;