Line data Source code
1 : /*----------------------------------------------------------------------------*/ 2 : /* CP2K: A general program to perform molecular dynamics simulations */ 3 : /* Copyright 2000-2023 CP2K developers group <https://cp2k.org> */ 4 : /* */ 5 : /* SPDX-License-Identifier: BSD-3-Clause */ 6 : /*----------------------------------------------------------------------------*/ 7 : 8 : #include "grid_sphere_cache.h" 9 : #include "grid_common.h" 10 : #include "grid_library.h" 11 : #include <assert.h> 12 : #include <math.h> 13 : #include <stdio.h> 14 : #include <stdlib.h> 15 : #include <string.h> 16 : 17 : /******************************************************************************* 18 : * \brief Compute the sphere bounds for a given single radius. 19 : * \author Ole Schuett 20 : ******************************************************************************/ 21 1534636 : static int single_sphere_bounds(const double disr_radius, const double dh[3][3], 22 : const double dh_inv[3][3], int *bounds) { 23 : 24 1534636 : int ibound = 0; 25 : 26 : // The cube contains an even number of grid points in each direction and 27 : // collocation is always performed on a pair of two opposing grid points. 28 : // Hence, the points with index 0 and 1 are both assigned distance zero via 29 : // the formular distance=(2*index-1)/2. 30 1534636 : const int kgmin = ceil(-1e-8 - disr_radius * dh_inv[2][2]); 31 1534636 : if (bounds != NULL) { 32 767318 : bounds[ibound] = kgmin; 33 : } 34 : ibound++; 35 15104632 : for (int kg = kgmin; kg <= 0; kg++) { 36 13569996 : const int kd = (2 * kg - 1) / 2; // distance from center in grid points 37 13569996 : const double kr = kd * dh[2][2]; // distance from center in a.u. 38 13569996 : const double kremain = disr_radius * disr_radius - kr * kr; 39 13569996 : const int jgmin = ceil(-1e-8 - sqrt(fmax(0.0, kremain)) * dh_inv[1][1]); 40 13569996 : if (bounds != NULL) { 41 6784998 : bounds[ibound] = jgmin; 42 : } 43 13569996 : ibound++; 44 141692618 : for (int jg = jgmin; jg <= 0; jg++) { 45 128122622 : const int jd = (2 * jg - 1) / 2; // distance from center in grid points 46 128122622 : const double jr = jd * dh[1][1]; // distance from center in a.u. 47 128122622 : const double jremain = kremain - jr * jr; 48 128122622 : const int igmin = ceil(-1e-8 - sqrt(fmax(0.0, jremain)) * dh_inv[0][0]); 49 128122622 : if (bounds != NULL) { 50 64061311 : bounds[ibound] = igmin; 51 : } 52 128122622 : ibound++; 53 : } 54 : } 55 1534636 : return ibound; // Number of bounds - needed to allocate array. 56 : } 57 : 58 : /******************************************************************************* 59 : * \brief Rebuild a cache entry for a given cell and max radius. 60 : * \author Ole Schuett 61 : ******************************************************************************/ 62 61199 : static void rebuild_cache_entry(const int max_imr, const double drmin, 63 : const double dh[3][3], 64 : const double dh_inv[3][3], 65 : grid_sphere_cache_entry *entry) { 66 61199 : if (entry->max_imr > 0) { 67 48097 : free(entry->offsets); 68 48097 : free(entry->storage); 69 : } 70 61199 : entry->max_imr = max_imr; 71 : 72 : // Compute required storage size. 73 61199 : entry->offsets = malloc(max_imr * sizeof(int)); 74 61199 : int nbounds_total = 0; 75 828517 : for (int imr = 1; imr <= max_imr; imr++) { 76 767318 : const double radius = imr * drmin; 77 767318 : const int nbounds = single_sphere_bounds(radius, dh, dh_inv, NULL); 78 767318 : entry->offsets[imr - 1] = nbounds_total; 79 767318 : nbounds_total += nbounds; 80 : } 81 : 82 : // Allocate and fill storage. 83 61199 : entry->storage = malloc(nbounds_total * sizeof(int)); 84 828517 : for (int imr = 1; imr <= max_imr; imr++) { 85 767318 : const double radius = imr * drmin; 86 767318 : const int offset = entry->offsets[imr - 1]; 87 767318 : single_sphere_bounds(radius, dh, dh_inv, &entry->storage[offset]); 88 : } 89 61199 : } 90 : 91 : /******************************************************************************* 92 : * \brief Lookup the sphere bound from cache and compute them as needed. 93 : * See grid_sphere_cache.h for details. 94 : * \author Ole Schuett 95 : ******************************************************************************/ 96 118958913 : void grid_sphere_cache_lookup(const double radius, const double dh[3][3], 97 : const double dh_inv[3][3], int **sphere_bounds, 98 : double *discr_radius) { 99 : 100 : // Prepare the cache. 101 118958913 : grid_sphere_cache *cache = grid_library_get_sphere_cache(); 102 : 103 : // Find or create cache entry for given grid. 104 118958913 : const double dr0 = dh[0][0], dr1 = dh[1][1], dr2 = dh[2][2]; 105 118958913 : grid_sphere_cache_entry *entry; 106 118958913 : bool found = false; 107 : 108 : // Fast path: check prev match. 109 118958913 : if (cache->prev_match < cache->size) { 110 118954584 : entry = &cache->entries[cache->prev_match]; 111 118954584 : if (entry->dr[0] == dr0 && entry->dr[1] == dr1 && entry->dr[2] == dr2) { 112 117428218 : found = true; 113 : } 114 : } 115 : 116 : // Full search. 117 118958913 : if (!found) { 118 4019201 : for (int i = 0; i < cache->size; i++) { 119 4006099 : entry = &cache->entries[i]; 120 4006099 : if (entry->dr[0] == dr0 && entry->dr[1] == dr1 && entry->dr[2] == dr2) { 121 1517593 : cache->prev_match = i; 122 1517593 : found = true; 123 1517593 : break; 124 : } 125 : } 126 : } 127 : 128 : // If no existing cache entry was found then create a new one. 129 118958913 : if (!found) { 130 13102 : cache->size++; 131 13102 : grid_sphere_cache_entry *old_entries = cache->entries; 132 13102 : const size_t entry_size = sizeof(grid_sphere_cache_entry); 133 13102 : cache->entries = malloc(cache->size * entry_size); 134 13102 : memcpy(cache->entries, old_entries, (cache->size - 1) * entry_size); 135 13102 : free(old_entries); 136 13102 : cache->prev_match = cache->size - 1; 137 13102 : entry = &cache->entries[cache->size - 1]; 138 : // Initialize new cache entry 139 13102 : entry->max_imr = 0; 140 13102 : entry->dr[0] = dr0; 141 13102 : entry->dr[1] = dr1; 142 13102 : entry->dr[2] = dr2; 143 13102 : entry->drmin = fmin(dr0, fmin(dr1, dr2)); 144 13102 : entry->drmin_inv = 1.0 / entry->drmin; 145 : } 146 : 147 : // Discretize the radius. 148 118958913 : const int imr = imax(1, (int)ceil(radius * entry->drmin_inv)); 149 118958913 : *discr_radius = entry->drmin * imr; 150 : 151 : // Rebuild cache entry if requested radius is too large. 152 118958913 : if (entry->max_imr < imr) { 153 61199 : rebuild_cache_entry(imr, entry->drmin, dh, dh_inv, entry); 154 : } 155 118958913 : const int offset = entry->offsets[imr - 1]; 156 118958913 : *sphere_bounds = &entry->storage[offset]; 157 118958913 : } 158 : 159 : /******************************************************************************* 160 : * \brief Free the memory of the sphere cache. 161 : * \author Ole Schuett 162 : ******************************************************************************/ 163 8266 : void grid_sphere_cache_free(grid_sphere_cache *cache) { 164 21368 : for (int i = 0; i < cache->size; i++) { 165 13102 : if (cache->entries[i].max_imr > 0) { 166 13102 : free(cache->entries[i].offsets); 167 13102 : free(cache->entries[i].storage); 168 : } 169 : } 170 8266 : free(cache->entries); 171 8266 : cache->size = 0; 172 8266 : } 173 : 174 : // EOF