Line data Source code
1 : /*----------------------------------------------------------------------------*/ 2 : /* CP2K: A general program to perform molecular dynamics simulations */ 3 : /* Copyright 2000-2024 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 1584888 : 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 1584888 : 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 1584888 : const int kgmin = ceil(-1e-8 - disr_radius * dh_inv[2][2]); 31 1584888 : if (bounds != NULL) { 32 792444 : bounds[ibound] = kgmin; 33 : } 34 : ibound++; 35 15584874 : for (int kg = kgmin; kg <= 0; kg++) { 36 13999986 : const int kd = (2 * kg - 1) / 2; // distance from center in grid points 37 13999986 : const double kr = kd * dh[2][2]; // distance from center in a.u. 38 13999986 : const double kremain = disr_radius * disr_radius - kr * kr; 39 13999986 : const int jgmin = ceil(-1e-8 - sqrt(fmax(0.0, kremain)) * dh_inv[1][1]); 40 13999986 : if (bounds != NULL) { 41 6999993 : bounds[ibound] = jgmin; 42 : } 43 13999986 : ibound++; 44 145902622 : for (int jg = jgmin; jg <= 0; jg++) { 45 131902636 : const int jd = (2 * jg - 1) / 2; // distance from center in grid points 46 131902636 : const double jr = jd * dh[1][1]; // distance from center in a.u. 47 131902636 : const double jremain = kremain - jr * jr; 48 131902636 : const int igmin = ceil(-1e-8 - sqrt(fmax(0.0, jremain)) * dh_inv[0][0]); 49 131902636 : if (bounds != NULL) { 50 65951318 : bounds[ibound] = igmin; 51 : } 52 131902636 : ibound++; 53 : } 54 : } 55 1584888 : 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 63191 : 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 63191 : if (entry->max_imr > 0) { 67 49562 : free(entry->offsets); 68 49562 : free(entry->storage); 69 : } 70 63191 : entry->max_imr = max_imr; 71 : 72 : // Compute required storage size. 73 63191 : entry->offsets = malloc(max_imr * sizeof(int)); 74 63191 : int nbounds_total = 0; 75 855635 : for (int imr = 1; imr <= max_imr; imr++) { 76 792444 : const double radius = imr * drmin; 77 792444 : const int nbounds = single_sphere_bounds(radius, dh, dh_inv, NULL); 78 792444 : entry->offsets[imr - 1] = nbounds_total; 79 792444 : nbounds_total += nbounds; 80 : } 81 : 82 : // Allocate and fill storage. 83 63191 : entry->storage = malloc(nbounds_total * sizeof(int)); 84 855635 : for (int imr = 1; imr <= max_imr; imr++) { 85 792444 : const double radius = imr * drmin; 86 792444 : const int offset = entry->offsets[imr - 1]; 87 792444 : single_sphere_bounds(radius, dh, dh_inv, &entry->storage[offset]); 88 : } 89 63191 : } 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 120216757 : 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 120216757 : grid_sphere_cache *cache = grid_library_get_sphere_cache(); 102 : 103 : // Find or create cache entry for given grid. 104 120216757 : const double dr0 = dh[0][0], dr1 = dh[1][1], dr2 = dh[2][2]; 105 120216757 : grid_sphere_cache_entry *entry = 0; 106 120216757 : bool found = false; 107 : 108 : // Fast path: check prev match. 109 120216757 : if (cache->prev_match < cache->size) { 110 120212272 : entry = &cache->entries[cache->prev_match]; 111 120212272 : if (entry->dr[0] == dr0 && entry->dr[1] == dr1 && entry->dr[2] == dr2) { 112 118658184 : found = true; 113 : } 114 : } 115 : 116 : // Full search. 117 120216757 : if (!found) { 118 4080149 : for (int i = 0; i < cache->size; i++) { 119 4066520 : entry = &cache->entries[i]; 120 4066520 : if (entry->dr[0] == dr0 && entry->dr[1] == dr1 && entry->dr[2] == dr2) { 121 1544944 : cache->prev_match = i; 122 1544944 : found = true; 123 1544944 : break; 124 : } 125 : } 126 : } 127 : 128 : // If no existing cache entry was found then create a new one. 129 120216757 : if (!found) { 130 13629 : cache->size++; 131 13629 : grid_sphere_cache_entry *old_entries = cache->entries; 132 13629 : const size_t entry_size = sizeof(grid_sphere_cache_entry); 133 13629 : cache->entries = malloc(cache->size * entry_size); 134 13629 : memcpy(cache->entries, old_entries, (cache->size - 1) * entry_size); 135 13629 : free(old_entries); 136 13629 : cache->prev_match = cache->size - 1; 137 13629 : entry = &cache->entries[cache->size - 1]; 138 : // Initialize new cache entry 139 13629 : entry->max_imr = 0; 140 13629 : entry->dr[0] = dr0; 141 13629 : entry->dr[1] = dr1; 142 13629 : entry->dr[2] = dr2; 143 13629 : entry->drmin = fmin(dr0, fmin(dr1, dr2)); 144 13629 : entry->drmin_inv = 1.0 / entry->drmin; 145 : } 146 : 147 : // Discretize the radius. 148 120216757 : const int imr = imax(1, (int)ceil(radius * entry->drmin_inv)); 149 120216757 : *discr_radius = entry->drmin * imr; 150 : 151 : // Rebuild cache entry if requested radius is too large. 152 120216757 : if (entry->max_imr < imr) { 153 63191 : rebuild_cache_entry(imr, entry->drmin, dh, dh_inv, entry); 154 : } 155 120216757 : const int offset = entry->offsets[imr - 1]; 156 120216757 : *sphere_bounds = &entry->storage[offset]; 157 120216757 : } 158 : 159 : /******************************************************************************* 160 : * \brief Free the memory of the sphere cache. 161 : * \author Ole Schuett 162 : ******************************************************************************/ 163 8430 : void grid_sphere_cache_free(grid_sphere_cache *cache) { 164 22059 : for (int i = 0; i < cache->size; i++) { 165 13629 : if (cache->entries[i].max_imr > 0) { 166 13629 : free(cache->entries[i].offsets); 167 13629 : free(cache->entries[i].storage); 168 : } 169 : } 170 8430 : free(cache->entries); 171 8430 : cache->size = 0; 172 8430 : } 173 : 174 : // EOF