1 static char help[] = "Test DMPlexGetLETKFLocalizationMatrix.\n\n"; 2 3 #include <petscdmplex.h> 4 5 int main(int argc, char **argv) 6 { 7 DM dm; 8 Mat H, Q; 9 PetscInt numobservations; 10 PetscInt dim = 1, n; 11 PetscInt faces[3]; 12 PetscReal lower[3] = {0.0, 0.0, 0.0}; 13 PetscReal upper[3] = {1.0, 1.0, 1.0}; 14 15 PetscFunctionBeginUser; 16 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 17 18 /* Get dimension and from options. We need the data here and Plex does not have access functions */ 19 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dm_plex_dim", &dim, NULL)); 20 n = 3; 21 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &n, NULL)); 22 n = 3; 23 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_lower", lower, &n, NULL)); 24 n = 3; 25 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_upper", upper, &n, NULL)); 26 27 /* Create the mesh using DMPlexCreateBoxMesh (could pass parameters) */ 28 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 29 PetscCall(DMSetType(dm, DMPLEX)); 30 PetscCall(DMSetFromOptions(dm)); 31 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 32 33 /* Verify dimension matches */ 34 PetscInt dmDim; 35 PetscCall(DMGetDimension(dm, &dmDim)); 36 PetscCheck(dmDim == dim, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "DM dimension %" PetscInt_FMT " does not match requested dimension %" PetscInt_FMT, dmDim, dim); 37 38 /* Set number of local observations to use: 3^dim */ 39 numobservations = 1; 40 for (PetscInt d = 0; d < dim && d < 2; d++) numobservations *= 3; 41 42 /* Get number of vertices */ 43 PetscInt vStart, vEnd, numVertices; 44 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 45 numVertices = vEnd - vStart; 46 47 /* Create a section for vertices (required for Global Point mapping) */ 48 PetscSection section; 49 PetscInt pStart, pEnd; 50 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 51 PetscCall(PetscSectionCreate(PETSC_COMM_WORLD, §ion)); 52 PetscCall(PetscSectionSetNumFields(section, 1)); 53 PetscCall(PetscSectionSetChart(section, pStart, pEnd)); 54 for (PetscInt v = vStart; v < vEnd; ++v) PetscCall(PetscSectionSetDof(section, v, 1)); 55 PetscCall(PetscSectionSetUp(section)); 56 PetscCall(DMSetLocalSection(dm, section)); 57 PetscCall(PetscSectionDestroy(§ion)); 58 59 /* Create global section */ 60 PetscSection globalSection; 61 PetscCall(DMGetGlobalSection(dm, &globalSection)); 62 63 /* Count observations (every other vertex in each dimension) */ 64 PetscInt numlocalobs = 0; 65 { 66 Vec coordinates; 67 PetscSection coordSection; 68 const PetscScalar *coordArray; 69 PetscInt offset; 70 71 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 72 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 73 PetscCall(VecGetArrayRead(coordinates, &coordArray)); 74 75 for (PetscInt v = vStart; v < vEnd; v++) { 76 PetscReal coords[3] = {0.0, 0.0, 0.0}; 77 PetscBool isObs = PETSC_TRUE; 78 79 PetscCall(PetscSectionGetOffset(coordSection, v, &offset)); 80 for (PetscInt d = 0; d < dim; d++) coords[d] = PetscRealPart(coordArray[offset + d]); 81 82 /* Check if this vertex is at an observation location (every other grid point) */ 83 for (PetscInt d = 0; d < dim; d++) { 84 PetscReal gridSpacing = (upper[d] - lower[d]) / faces[d]; 85 PetscInt gridIdx = (PetscInt)((coords[d] - lower[d]) / gridSpacing + 0.5); 86 if (gridIdx % 2 != 0) { 87 isObs = PETSC_FALSE; 88 break; 89 } 90 } 91 if (isObs) numlocalobs++; 92 } 93 PetscCall(VecRestoreArrayRead(coordinates, &coordArray)); 94 } 95 96 /* Create H matrix */ 97 PetscCall(MatCreate(PETSC_COMM_WORLD, &H)); 98 PetscCall(MatSetSizes(H, numlocalobs, PETSC_DECIDE, PETSC_DECIDE, numVertices)); 99 PetscCall(MatSetType(H, MATAIJ)); 100 PetscCall(MatSeqAIJSetPreallocation(H, 1, NULL)); 101 PetscCall(MatMPIAIJSetPreallocation(H, 1, NULL, 0, NULL)); 102 PetscCall(PetscObjectSetName((PetscObject)H, "H_observation_operator")); 103 104 /* Fill H matrix */ 105 { 106 Vec coordinates; 107 PetscSection coordSection; 108 const PetscScalar *coordArray; 109 PetscInt obsIdx = 0; 110 PetscInt offset; 111 112 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 113 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 114 PetscCall(VecGetArrayRead(coordinates, &coordArray)); 115 116 for (PetscInt v = vStart; v < vEnd; v++) { 117 PetscReal coords[3] = {0.0, 0.0, 0.0}; 118 PetscBool isObs = PETSC_TRUE; 119 120 PetscCall(PetscSectionGetOffset(coordSection, v, &offset)); 121 for (PetscInt d = 0; d < dim; d++) coords[d] = PetscRealPart(coordArray[offset + d]); 122 123 /* Check if this vertex is at an observation location (every other grid point) */ 124 for (PetscInt d = 0; d < dim; d++) { 125 PetscReal gridSpacing = (upper[d] - lower[d]) / faces[d]; 126 PetscInt gridIdx = (PetscInt)((coords[d] - lower[d]) / gridSpacing + 0.5); 127 if (gridIdx % 2 != 0) { 128 isObs = PETSC_FALSE; 129 break; 130 } 131 } 132 133 if (isObs) { 134 PetscCall(MatSetValue(H, obsIdx, v - vStart, 1.0, INSERT_VALUES)); 135 obsIdx++; 136 } 137 } 138 PetscCall(VecRestoreArrayRead(coordinates, &coordArray)); 139 } 140 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 141 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 142 143 /* View H */ 144 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Observation Operator H:\n")); 145 PetscCall(MatView(H, PETSC_VIEWER_STDOUT_WORLD)); 146 147 /* Perturb interior vertex coordinates */ 148 { 149 Vec coordinates; 150 PetscSection coordSection; 151 PetscScalar *coordArray; 152 unsigned long seed = 123456789; 153 154 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 155 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 156 PetscCall(VecGetArray(coordinates, &coordArray)); 157 158 for (PetscInt v = vStart; v < vEnd; v++) { 159 PetscInt offset; 160 161 PetscCall(PetscSectionGetOffset(coordSection, v, &offset)); 162 163 for (PetscInt d = 0; d < dim; d++) { 164 PetscReal noise, gridSpacing = (upper[d] - lower[d]) / faces[d]; 165 166 seed = (1103515245 * seed + 12345) % 2147483648; 167 noise = (PetscReal)seed / 2147483648.0; 168 coordArray[offset + d] += (noise - 0.5) * 0.05 * gridSpacing; 169 } 170 } 171 PetscCall(VecRestoreArray(coordinates, &coordArray)); 172 } 173 174 /* Call the function */ 175 PetscCall(DMPlexGetLETKFLocalizationMatrix(dm, numobservations, numlocalobs, H, &Q)); 176 PetscCall(PetscObjectSetName((PetscObject)Q, "Q_localization")); 177 178 /* View Q */ 179 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Localization Matrix Q:\n")); 180 PetscCall(MatView(Q, PETSC_VIEWER_STDOUT_WORLD)); 181 182 /* Cleanup */ 183 PetscCall(MatDestroy(&H)); 184 PetscCall(MatDestroy(&Q)); 185 PetscCall(DMDestroy(&dm)); 186 PetscCall(PetscFinalize()); 187 return 0; 188 } 189 190 /*TEST 191 192 test: 193 requires: kokkos 194 suffix: 1 195 diff_args: -j 196 args: -dm_plex_dim 1 -dm_plex_box_faces 16 -dm_plex_simplex 0 197 198 test: 199 requires: kokkos 200 suffix: 2 201 diff_args: -j 202 args: -dm_plex_dim 2 -dm_plex_box_faces 8,8 -dm_plex_simplex 0 203 204 test: 205 requires: kokkos 206 suffix: 3 207 diff_args: -j 208 args: -dm_plex_dim 3 -dm_plex_box_faces 5,5,5 -dm_plex_simplex 0 209 210 TEST*/ 211