1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
|
/*
* Pseudospectral solver for 2nd order 2D linear PDE systems
* Copyright (C) 2014-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/>.
*/
#ifndef TEUKOLSKY_PSSOLVE_H
#define TEUKOLSKY_PSSOLVE_H
/**
* The problem being solved is a sequence of N linear partial differential
* equations
*
* / \
* ∑ | ∑ C_{ab}^{ik} ∂_a ∂_b u_k + ∑ C_{a}^{ik} ∂_a u_k + C^{ik} u_k | = S^i
* k \ a,b a /
*
* where
* * i numbers the equations and runs from 0 to N-1
* * k numbers the unknown functions and also runs from 0 to N-1
* * a and b identify spatial directions and run from 0 to 1
* * u_k = u_k(x_a) is the k-th unknown function
* * C_{ab}^{ik}, C_{a}^{ik} and C^{ik} are the coefficients in front of
* the corresponding derivative of k-th unknown function in the i-th
* equation
* * S^i is the right-hand side of the i-th equation
* C_{*}^{ik} and S^i are all (known) functions of space and define the
* equations to be solved.
*/
#include "config.h"
#if HAVE_OPENCL
#include <cl.h>
#else
typedef void* cl_context;
typedef void* cl_command_queue;
#endif
#include <stdint.h>
#include <threadpool.h>
#include "basis.h"
#include "log.h"
enum PSSolveDiffOrder {
PSSOLVE_DIFF_ORDER_00,
PSSOLVE_DIFF_ORDER_10,
PSSOLVE_DIFF_ORDER_01,
PSSOLVE_DIFF_ORDER_11,
PSSOLVE_DIFF_ORDER_20,
PSSOLVE_DIFF_ORDER_02,
PSSOLVE_DIFF_ORDER_NB,
};
typedef struct PSSolvePriv PSSolvePriv;
typedef struct PSSolveContext {
/**
* Solver private data, not to be touched by the caller.
*/
PSSolvePriv *priv;
/**
* The logging context.
* Set by the caller before tdi_pssolve_context_init().
*/
TDLogger logger;
/**
* Number of equations/unknown functions in the set.
* Set by tdi_pssolve_context_alloc().
*/
unsigned int nb_equations;
/**
* The basis sets.
*
* basis[i][j] is the basis set used for i-th variable in j-th direction.
*
* The array is allocated by tdi_pssolve_context_alloc(), must be filled by
* by the caller before tdi_pssolve_context_init().
*/
const BasisSetContext *(*basis)[2];
/**
* Order of the solver.
*
* solve_order[i][j] is the order of the solver (i.e. the number of the
* basis functions used) for i-th variable in j-th direction.
*
* Allocated by tdi_pssolve_context_alloc(), must be filled by the caller
* before tdi_pssolve_context_init().
*/
unsigned int (*solve_order)[2];
/**
* Locations of the collocation points. The equation coefficients passed to
* tdi_pssolve_solve() should be evaluated at those grid positions.
*
* colloc_grid[i][j] is an array of length solve_order[i][j] and contains
* the collocation points for the i-th equation in the j-th direction.
*
* Set by the solver after tdi_pssolve_context_init().
*/
const double * const (*colloc_grid)[2];
/**
* The thread pool used for multithreaded execution. May be set by the
* caller before tdi_pssolve_context_init(), otherwise a single thread will
* be used.
*/
TPContext *tp;
uint64_t lu_solves_count;
uint64_t lu_solves_time;
uint64_t cg_solve_count;
uint64_t cg_iter_count;
uint64_t cg_time_total;
uint64_t construct_matrix_count;
uint64_t construct_matrix_time;
} PSSolveContext;
typedef struct PSEqCoeffs {
const double * const (*func_coeffs)[PSSOLVE_DIFF_ORDER_NB];
} PSEqCoeffs;
/**
* Allocate a new solver.
*
* @param ctx The newly allocated solver context will be written here.
* @param nb_equations number of equations to solve (equal to the number of
* unknown functions to solve for)
*
* @return 0 on success a negative error code on failure.
*/
int tdi_pssolve_context_alloc(PSSolveContext **ctx, unsigned int nb_equations);
/**
* Initialize the solver for use after all the context options have been set.
* This function must be called exactly once before any calls to
* tdi_pssolve_solve().
*
* @return 0 on success, a negative error code on failure.
*/
int tdi_pssolve_context_init(PSSolveContext *ctx);
/**
* Free the solver and write NULL to the supplied pointer.
*/
void tdi_pssolve_context_free(PSSolveContext **ctx);
/**
* Solve a PDE. This function may be called multiple times in succession to
* solve multiple related PDEs (it will be efficient if the equation
* coefficients do not change too much).
*
* @param ctx the solver context
* @param eq_coeffs the equation coefficients at the collocation points.
* eq_coeffs[i].func_coeffs[j][k] is the array of coefficients
* for the k-th derivative (as per enum PSSolveDiffOrder) of
* the j-th unknown function in the i-th equation.
* @param rhs the right-hand side of the equation at the collocation points.
* @param coeffs the spectral coefficients of the solution will be written here.
*
* @return 0 on success, a negative error code on failure. The contents of
* coeffs are undefined on failure.
*/
int tdi_pssolve_solve(PSSolveContext *ctx,
const PSEqCoeffs *eq_coeffs,
const double *rhs, double *coeffs);
int tdi_pssolve_diff_order(enum PSSolveDiffOrder order, unsigned int dir);
#endif /* TEUKOLSKY_PSSOLVE_H */
|
