diff --git a/src/hf/bsplines_and_grid.c b/src/hf/bsplines_and_grid.c index 318df67..ab8f2fe 100644 --- a/src/hf/bsplines_and_grid.c +++ b/src/hf/bsplines_and_grid.c @@ -1,16 +1,27 @@ -global BSplineCtx g_bspline_ctx = {0}; -global Grid g_grid = {0}; + +global BSplineCtx *g_bspline_ctx = 0; +global Grid *g_grid = 0; global U32 g_debug_bspline_matrix = 0; +function void +bspline_ctx_assign(BSplineCtx *ctx) { + g_bspline_ctx = ctx; +} + +function void +grid_assign(Grid *grid) { + g_grid = grid; +} + function F64 bspline_recursion(F64 x, U32 k, U32 i) { - F64 *t = g_bspline_ctx.knotpoints; + F64 *t = g_bspline_ctx->knotpoints; F64 tolerance = 1e-14; if(k == 1) { - if(i == g_bspline_ctx.num_bsplines-1 && fabs(x -g_grid.end) < tolerance ) + if(i == g_bspline_ctx->num_bsplines-1 && fabs(x - g_grid->end) < tolerance ) { // TODO(anton): // This is like a hack to get the last bspline to be 1 at the last point. @@ -18,7 +29,7 @@ bspline_recursion(F64 x, U32 k, U32 i) // to unity at the last point, actually. I have to check this. return 1.0; } - else if(i < g_bspline_ctx.num_bsplines && (x >= t[i] && x < t[i+1])) + else if(i < g_bspline_ctx->num_bsplines && (x >= t[i] && x < t[i+1])) { return 1.0; } else { @@ -45,7 +56,7 @@ bspline_recursion(F64 x, U32 k, U32 i) function F64 compute_bspline_F64(F64 x_coord, U32 index) { - U32 k = g_bspline_ctx.order; + U32 k = g_bspline_ctx->order; F64 out = bspline_recursion(x_coord, k, index); return out; } @@ -54,9 +65,9 @@ compute_bspline_F64(F64 x_coord, U32 index) function F64 compute_dBspline_F64(F64 x_coord, U32 index) { - U32 k = g_bspline_ctx.order; + U32 k = g_bspline_ctx->order; F64 prefac = (F64)(k - 1); - F64 *t = g_bspline_ctx.knotpoints; + F64 *t = g_bspline_ctx->knotpoints; F64 term1_enum = bspline_recursion(x_coord, k-1, index); F64 term2_enum = bspline_recursion(x_coord, k-1, index+1); F64 term1_denom = t[index+k-1] - t[index]; @@ -70,17 +81,17 @@ compute_dBspline_F64(F64 x_coord, U32 index) function void set_up_grid(Arena *arena) { - g_grid.start = GRID_START_POINT; - g_grid.end = GRID_END_POINT; - g_grid.num_steps = GRID_NUM_STEPS; + g_grid->start = GRID_START_POINT; + g_grid->end = GRID_END_POINT; + g_grid->num_steps = GRID_NUM_STEPS; - g_grid.points = PushArray(arena, F64, g_grid.num_steps); - F64 step_size = (g_grid.end-g_grid.start)/(F64)g_grid.num_steps; - g_grid.points[0] = g_grid.start; - g_grid.points[g_grid.num_steps-1] = g_grid.end; - for(U32 i = 1; i < g_grid.num_steps-1; i++) + g_grid->points = PushArray(arena, F64, g_grid->num_steps); + F64 step_size = (g_grid->end-g_grid->start)/(F64)g_grid->num_steps; + g_grid->points[0] = g_grid->start; + g_grid->points[g_grid->num_steps-1] = g_grid->end; + for(U32 i = 1; i < g_grid->num_steps-1; i++) { - g_grid.points[i] = g_grid.points[i-1] + step_size; + g_grid->points[i] = g_grid->points[i-1] + step_size; } } @@ -95,13 +106,13 @@ get_bspline_index_size(U32 size1, U32 i, U32 j) function inline U32 get_bspline_index(U32 i, U32 j) { - return g_bspline_ctx.num_knotpoints * j + i; + return g_bspline_ctx->num_knotpoints * j + i; } function inline U32 get_bspline_grid_index(U32 i, U32 j) { - return g_grid.num_steps * j + i; + return g_grid->num_steps * j + i; } @@ -111,29 +122,29 @@ set_up_bspline_context(Arena* arena) // Create knotpoint sequence. U32 k = 4; U32 bspl_N = 40; - g_bspline_ctx.order = k; - g_bspline_ctx.num_knotpoints = bspl_N; - g_bspline_ctx.num_bsplines = bspl_N-k; - g_bspline_ctx.num_phys_points = bspl_N-(2*k)+2; // Remove k points at each end, and then add back the first and last points. - g_bspline_ctx.arena = arena; - g_bspline_ctx.knotpoints = PushArray(arena, F64, g_bspline_ctx.num_knotpoints); + g_bspline_ctx->order = k; + g_bspline_ctx->num_knotpoints = bspl_N; + g_bspline_ctx->num_bsplines = bspl_N-k; + g_bspline_ctx->num_phys_points = bspl_N-(2*k)+2; // Remove k points at each end, and then add back the first and last points. + g_bspline_ctx->arena = arena; + g_bspline_ctx->knotpoints = PushArray(arena, F64, g_bspline_ctx->num_knotpoints); // Set up physical points; - F64 delta = (g_grid.end-g_grid.start)/(g_bspline_ctx.num_phys_points-1); + F64 delta = (g_grid->end-g_grid->start)/(g_bspline_ctx->num_phys_points-1); // Set ghost points including first physical - U32 phys_point_last_index = g_bspline_ctx.num_phys_points + k-1; + U32 phys_point_last_index = g_bspline_ctx->num_phys_points + k-1; for(U32 i = 0; i < k; i++) { - g_bspline_ctx.knotpoints[i] = g_grid.start; + g_bspline_ctx->knotpoints[i] = g_grid->start; } for(U32 i = k; i < phys_point_last_index; i++) { - g_bspline_ctx.knotpoints[i] = g_bspline_ctx.knotpoints[i-1] + delta; + g_bspline_ctx->knotpoints[i] = g_bspline_ctx->knotpoints[i-1] + delta; } // Set the last points - F64 last_physical = g_grid.end; - for(U32 i = phys_point_last_index; i < g_bspline_ctx.num_knotpoints; i++) + F64 last_physical = g_grid->end; + for(U32 i = phys_point_last_index; i < g_bspline_ctx->num_knotpoints; i++) { - g_bspline_ctx.knotpoints[i] = last_physical; + g_bspline_ctx->knotpoints[i] = last_physical; } } @@ -185,11 +196,11 @@ function void set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena) { - U64 num_bsplines = g_bspline_ctx.num_bsplines; - U64 k = g_bspline_ctx.order; - F64 *t = g_bspline_ctx.knotpoints; - U32 num_grid_points = g_grid.num_steps; - U32 num_knotpoints = g_bspline_ctx.num_knotpoints; + U64 num_bsplines = g_bspline_ctx->num_bsplines; + U64 k = g_bspline_ctx->order; + F64 *t = g_bspline_ctx->knotpoints; + U32 num_grid_points = g_grid->num_steps; + U32 num_knotpoints = g_bspline_ctx->num_knotpoints; // For sanity check we make the first 4 bsplines by hand. { @@ -205,7 +216,7 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena) F64 *dBspl9 = PushArray(arena, F64, num_grid_points); for(U32 i = 0; i < num_grid_points; i++) { - F64 x = g_grid.points[i]; + F64 x = g_grid->points[i]; bspl0[i] = compute_bspline_F64(x, 0); bspl1[i] = compute_bspline_F64(x, 1); bspl2[i] = compute_bspline_F64(x, 2); @@ -218,7 +229,7 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena) dBspl9[i] = compute_dBspline_F64(x, 9); } - F64 test = compute_bspline_F64(g_grid.points[num_grid_points-1], 9); + F64 test = compute_bspline_F64(g_grid->points[num_grid_points-1], 9); write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\bspline0.dat"), bspl0, num_grid_points, "%13.6e\n"); write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\bspline1.dat"), bspl1, num_grid_points, "%13.6e\n"); @@ -235,18 +246,18 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena) { ArenaTemp scratch = scratch_get(0,0); - g_bspline_ctx.bsplines_grid = PushArray(arena, F64, num_grid_points*num_bsplines); - g_bspline_ctx.dBsplines_grid = PushArray(arena, F64, num_grid_points*num_bsplines); + g_bspline_ctx->bsplines_grid = PushArray(arena, F64, num_grid_points*num_bsplines); + g_bspline_ctx->dBsplines_grid = PushArray(arena, F64, num_grid_points*num_bsplines); for(U32 j = 0; j < num_bsplines; j++) { for(U32 i = 0; i < num_grid_points; i++) { U32 index = get_bspline_grid_index(i, j); - F64 bspline_value = compute_bspline_F64(g_grid.points[i], j); - F64 dBspline_value = compute_dBspline_F64(g_grid.points[i], j); - g_bspline_ctx.bsplines_grid[index] = bspline_value; - g_bspline_ctx.dBsplines_grid[index] = dBspline_value; + F64 bspline_value = compute_bspline_F64(g_grid->points[i], j); + F64 dBspline_value = compute_dBspline_F64(g_grid->points[i], j); + g_bspline_ctx->bsplines_grid[index] = bspline_value; + g_bspline_ctx->dBsplines_grid[index] = dBspline_value; if(g_debug_bspline_matrix && j == 0) { String8 out = str8_pushf(scratch.arena, "%i %i \t %13.6e\n", i, index, bspline_value); @@ -258,35 +269,35 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena) scratch_release(scratch); } - g_bspline_ctx.bsplines = PushArray(arena, F64, num_knotpoints*num_bsplines); - g_bspline_ctx.dBsplines = PushArray(arena, F64, num_knotpoints*num_bsplines); + g_bspline_ctx->bsplines = PushArray(arena, F64, num_knotpoints*num_bsplines); + g_bspline_ctx->dBsplines = PushArray(arena, F64, num_knotpoints*num_bsplines); for(U32 i = 0; i < num_knotpoints; i++) { for(U32 j = 0; j < num_bsplines; j++) { U32 index = get_bspline_index(i, j); - g_bspline_ctx.bsplines[index] = compute_bspline_F64(t[i], j); - g_bspline_ctx.dBsplines[index] = compute_dBspline_F64(t[i], j); + g_bspline_ctx->bsplines[index] = compute_bspline_F64(t[i], j); + g_bspline_ctx->dBsplines[index] = compute_dBspline_F64(t[i], j); } } - write_bspline_matrix_F64(g_bspline_ctx.bsplines_grid, + write_bspline_matrix_F64(g_bspline_ctx->bsplines_grid, num_grid_points, num_bsplines, str8_lit(bspline_grid_array_file_path)); - write_bspline_matrix_F64(g_bspline_ctx.dBsplines_grid, + write_bspline_matrix_F64(g_bspline_ctx->dBsplines_grid, num_grid_points, num_bsplines, str8_lit(dBspline_grid_array_file_path)); - write_bspline_matrix_F64(g_bspline_ctx.bsplines, + write_bspline_matrix_F64(g_bspline_ctx->bsplines, num_knotpoints, num_bsplines, str8_lit(bspline_knots_array_file_path)); - write_bspline_matrix_F64(g_bspline_ctx.dBsplines, + write_bspline_matrix_F64(g_bspline_ctx->dBsplines, num_knotpoints, num_bsplines, str8_lit(dBspline_knots_array_file_path)); diff --git a/src/hf/bsplines_and_grid.h b/src/hf/bsplines_and_grid.h index 02f4314..13863f0 100644 --- a/src/hf/bsplines_and_grid.h +++ b/src/hf/bsplines_and_grid.h @@ -41,5 +41,6 @@ function inline U32 get_bspline_grid_index(U32 i, U32 j); function void set_up_bspline_context(Arena* arena); function void write_bspline_matrix_F64(F64 *bsplines, U32 size1, U32 size2, String8 filename_path); function void set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena); - +function void bspline_ctx_assign(BSplineCtx *ctx); +function void grid_assign(Grid *grid); #endif /* BSPLINES_AND_GRID_H */ diff --git a/src/hf/hf_base.c b/src/hf/hf_base.c index 2bdd342..31c35f2 100644 --- a/src/hf/hf_base.c +++ b/src/hf/hf_base.c @@ -1,6 +1,5 @@ - //~ Globals -global GaussLegendre g_gauss_legendre = {0}; +GaussLegendre g_gauss_legendre = {0}; function inline U32 mat_get_col_major_idx(U32 i, U32 j, U32 size1) { diff --git a/src/main.c b/src/main.c index 824c05b..ede292c 100644 --- a/src/main.c +++ b/src/main.c @@ -29,30 +29,75 @@ typedef struct Eigensolution_F64 Eigensolution_F64; struct Eigensolution_F64 { U32 l; - F64 *eigenvalues_re; - F64 *eigenvalues_im; - Mat_F64 right_eigenvectors; - Mat_F64 left_eigenvectors; + F64 *eigvals_re; + F64 *eigvals_im; + Mat_F64 right_eigvecs; + Mat_F64 left_eigvecs; +}; + +typedef enum ElectronConfig ElectronConfig; +enum ElectronConfig { + ECFG_1s, + ECFG_2s, + ECFG_2p, + ECFG_3s, + ECFG_3p, + ECFG_3d, + ECFG_4s, + ECFG_4p, + ECFG_4d, + ECFG_4f, + ECFG_NUM_CONFIGS +}; + +typedef enum AngularMomenta AngularMomenta; +enum AngularMomenta { + ANGMOM_s, + ANGMOM_p, + ANGMOM_d, + ANGMOM_f, + ANGMOM_NUM_MOMENTA }; -#define NUM_ANGULAR_MOMENTA 3 typedef struct Atom Atom; struct Atom { - Arena *arena; + String8 name; U32 Z; - Eigensolution_F64 eigensolutions[NUM_ANGULAR_MOMENTA]; + U32 occupancy[ECFG_NUM_CONFIGS]; +}; + + +// We use a "fat struct" approach where everything just exists here +// in a single struct. +#define MAX_NUM_ANGULAR_MOMENTA 3 +typedef struct Problem Problem; +struct Problem { + Arena *arena; // Just use a single arena to start with + Grid grid; + BSplineCtx bspline_ctx; + Atom atom; + F64 angular_momentum_l[ANGMOM_NUM_MOMENTA]; + Eigensolution_F64 eigsols[MAX_NUM_ANGULAR_MOMENTA]; + U32 num_eigsols; + Mat_F64 H; + Mat_F64 H_l; }; ////// //~ -global Arena *g_base_arena = 0; -global Arena *g_filename_arena = 0; -global Atom g_atom = {0}; -global F64 g_angular_momenta[NUM_ANGULAR_MOMENTA] = {0.0, 1.0, 2.0}; - -////// -//~ +function Problem problem_create() { + Problem out = {0}; + out.arena = m_make_arena(); + out.atom.name = str8_lit("Hydrogen"); + out.atom.Z = 1; + out.atom.occupancy[ECFG_1s] = 1; + out.angular_momentum_l[ANGMOM_s] = 0.0; + out.angular_momentum_l[ANGMOM_p] = 1.0; + out.angular_momentum_l[ANGMOM_d] = 2.0; + out.angular_momentum_l[ANGMOM_f] = 3.0; + return out; +} /* Auxiliary routine: printing a matrix */ function void print_eigenvalues(S32 l, S32 n, F64 *wr, F64 *wi) { @@ -77,16 +122,16 @@ function void print_eigenvalues(S32 l, S32 n, F64 *wr, F64 *wi) { } function void -set_up_first_matrices(Atom *atom, Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) { +set_up_first_matrices(Problem *problem, Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) { // We work in units hbar = 1, bohr radius a0 = 1, electron mass m_e = 1, and // charge e = 1, and 1/(4piepsilon_0) = 1. Set up Hamiltonian: H = // -0.5*d^2/dr^2 + l(l+1)/(2r^2) - Z/r { ArenaTemp scratch = scratch_get(0, 0); - - F64 *t = g_bspline_ctx.knotpoints; - F64 Z = (F64)atom->Z; - U32 k = g_bspline_ctx.order; + BSplineCtx *bspl_ctx = &problem->bspline_ctx; + F64 *t = bspl_ctx->knotpoints; + F64 Z = (F64)problem->atom.Z; + U32 k = bspl_ctx->order; // Skipping first bspline for (U32 i = 0; i < H->size1; i++) { @@ -159,7 +204,7 @@ set_up_first_matrices(Atom *atom, Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) { } LOG(str8_pushf(scratch.arena, "H.size1=N-k-2=%i, last bspline index=%i \n", - H->size1, g_bspline_ctx.num_bsplines - 1) + H->size1, bspl_ctx->num_bsplines - 1) .str); scratch_release(scratch); //print_mat_F64(H); @@ -174,9 +219,9 @@ function void compute_wf_norm_F64(F64 *coeffs, U64 coeff_size, U64 n, U64 l) { // Gauss legendre integration // F64 norm = 0.0; - for (U64 i = 0; i < g_grid.num_steps - 1; i++) { - F64 a = g_grid.points[i]; - F64 b = g_grid.points[i + 1]; + for (U64 i = 0; i < g_grid->num_steps - 1; i++) { + F64 a = g_grid->points[i]; + F64 b = g_grid->points[i + 1]; F64 prefac = 0.5 * (b - a); // Only integrate non-zero intervals @@ -209,46 +254,45 @@ function void EntryPoint(void) { OS_InitReceipt os_receipt = OS_init(); OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt); - - g_filename_arena = m_make_arena_reserve(Megabytes(2)); - g_base_arena = m_make_arena(); - g_atom.Z = 1.0; // Hydrogen - g_atom.arena = m_make_arena(); + Problem problem = problem_create(); + LOG(str8_pushf(problem.arena, "Created Problem-struct for %s \n", problem.atom.name).str); - set_up_gauss_legendre_points(g_base_arena); + set_up_gauss_legendre_points(problem.arena); //- Set up grid and write to file. - set_up_grid(g_base_arena); + grid_assign(&problem.grid); + set_up_grid(problem.arena); write_array_binary_F64(str8_lit(grid_file_path_bin), - g_grid.points, - g_grid.num_steps); - write_array_F64(str8_lit(grid_file_path), g_grid.points, g_grid.num_steps, + problem.grid.points, + problem.grid.num_steps); + write_array_F64(str8_lit(grid_file_path), problem.grid.points,problem.grid.num_steps, "%13.6e\n"); //- The BSpline context is the knotpoints and the BSpline order etc. - set_up_bspline_context(g_base_arena); - write_array_F64(str8_lit(knotpoints_file_path), g_bspline_ctx.knotpoints, - g_bspline_ctx.num_knotpoints, "%13.6e\n"); + bspline_ctx_assign(&problem.bspline_ctx); + set_up_bspline_context(problem.arena); + write_array_F64(str8_lit(knotpoints_file_path), problem.bspline_ctx.knotpoints, + problem.bspline_ctx.num_knotpoints, "%13.6e\n"); //- Then we generate the BSplines and save them off for reference and // debugging. - set_up_bsplines_at_points_and_write_matrix_F64(g_base_arena); + set_up_bsplines_at_points_and_write_matrix_F64(problem.arena); - U32 N = g_bspline_ctx.num_knotpoints; - U32 k = g_bspline_ctx.order; + U32 N = problem.bspline_ctx.num_knotpoints; + U32 k = problem.bspline_ctx.order; U32 mat_size1 = N - k - 2; U32 mat_size2 = mat_size1; - Mat_F64 H_base = mat_F64(g_base_arena, mat_size1, mat_size2); - Mat_F64 H_l_base = mat_F64(g_base_arena, mat_size1, mat_size2); - Mat_F64 H_l = mat_F64(g_base_arena, mat_size1, mat_size2); - Mat_F64 H = mat_F64(g_base_arena, mat_size1, mat_size2); + Mat_F64 H_base = mat_F64(problem.arena, mat_size1, mat_size2); + Mat_F64 H_l_base = mat_F64(problem.arena, mat_size1, mat_size2); + Mat_F64 H_l = mat_F64(problem.arena, mat_size1, mat_size2); + Mat_F64 H = mat_F64(problem.arena, mat_size1, mat_size2); // This will be the inverse of B, but to start with we construct B. - Mat_F64 B_inv = mat_F64(g_base_arena, mat_size1, mat_size2); + Mat_F64 B_inv = mat_F64(problem.arena, mat_size1, mat_size2); // A is the actual matrix for each eigenvalue problem. - Mat_F64 A = mat_F64(g_base_arena, H.size1, H.size2); - set_up_first_matrices(&g_atom, &H_base, &H_l_base, &B_inv); + Mat_F64 A = mat_F64(problem.arena, H.size1, H.size2); + set_up_first_matrices(&problem, &H_base, &H_l_base, &B_inv); // Our problem is Hc = EBc, but we want to solve B^-1Hc = Ec, // so we invert the B matrix and compute the product A = B^-1H before calling @@ -257,11 +301,11 @@ function void EntryPoint(void) { // This arena is used to push results from f. ex eigenvalue computations. // For each angular momentum - for (U32 ang_mom_idx = 0; ang_mom_idx < NUM_ANGULAR_MOMENTA; ang_mom_idx++) { + for (U32 ang_mom_idx = 0; ang_mom_idx < MAX_NUM_ANGULAR_MOMENTA; ang_mom_idx++) { ArenaTemp scratch = scratch_get(0, 0); mat_F64_copy_to_dst(&H, &H_base); - F64 l = g_angular_momenta[ang_mom_idx]; - Eigensolution_F64 *eigsol = &g_atom.eigensolutions[ang_mom_idx]; + F64 l = problem.angular_momentum_l[ang_mom_idx]; + Eigensolution_F64 *eigsol = &problem.eigsols[ang_mom_idx]; eigsol->l = (U32)l; if (l > 1e-16) { @@ -295,14 +339,14 @@ function void EntryPoint(void) { F64 wkopt; F64 *work; - eigsol->eigenvalues_re = PushArray(g_atom.arena, F64, size1); - F64 *wr = eigsol->eigenvalues_re; - eigsol->eigenvalues_im = PushArray(g_atom.arena, F64, size1); - F64 *wi = eigsol->eigenvalues_im; - eigsol->left_eigenvectors = mat_F64(g_atom.arena, ldvl, size1); - F64 *vl = eigsol->left_eigenvectors.data; - eigsol->right_eigenvectors = mat_F64(g_atom.arena, size1, ldvr); - F64 *vr = eigsol->right_eigenvectors.data; + eigsol->eigvals_re= PushArray(problem.arena, F64, size1); + F64 *wr = eigsol->eigvals_re; + eigsol->eigvals_im = PushArray(problem.arena, F64, size1); + F64 *wi = eigsol->eigvals_im; + eigsol->left_eigvecs = mat_F64(problem.arena, ldvl, size1); + F64 *vl = eigsol->left_eigvecs.data; + eigsol->right_eigvecs = mat_F64(problem.arena, size1, ldvr); + F64 *vr = eigsol->right_eigvecs.data; lwork = -1; F64 *a = A.data; @@ -339,9 +383,9 @@ function void EntryPoint(void) { // size1, n, ang_mom_idx); U64 eigvec_idx = mat_get_col_major_idx(0, energy_index, size1); - F64 *eigvecs = &eigsol->right_eigenvectors.data[eigvec_idx]; + F64 *eigvecs = &eigsol->right_eigvecs.data[eigvec_idx]; write_array_F64( - get_eigenvector_filename(g_filename_arena, n, ang_mom_idx), + get_eigenvector_filename(scratch.arena, n, ang_mom_idx), eigvecs, size1, "%13.6e\n");