PSCF v1.3
rpg/field/FieldIo.tpp
1#ifndef RPG_FIELD_IO_TPP
2#define RPG_FIELD_IO_TPP
3
4/*
5* PSCF - Polymer Self-Consistent Field
6*
7* Copyright 2015 - 2025, The Regents of the University of Minnesota
8* Distributed under the terms of the GNU General Public License.
9*/
10
11#include "FieldIo.h"
12
13#include <rpg/field/HostDArrayComplex.h>
14
15#include <pscf/math/complex.h>
16
17#include <prdc/field/FieldIoReal.tpp> // base class implementation
18#include <prdc/field/fieldIoUtil.h>
19#include <prdc/field/fieldArrayUtil.h>
20#include <prdc/crystal/Basis.h>
21#include <prdc/crystal/UnitCell.h>
22#include <prdc/cuda/RFieldComparison.h>
23#include <prdc/cuda/types.h>
24#include <prdc/cuda/complex.h>
25#include <prdc/cuda/resources.h>
26
27#include <pscf/mesh/Mesh.h>
28#include <pscf/math/IntVec.h>
29#include <pscf/cuda/HostDArray.h>
30
31
32namespace Pscf {
33namespace Rpg {
34
35 using namespace Util;
36 using namespace Pscf::Prdc;
37 using namespace Pscf::Prdc::Cuda;
38
39 // Field Io in r-grid format
40
41 /*
42 * Read an array of fields in r-grid format.
43 */
44 template <int D>
46 std::istream &in,
47 DArray<RField<D> >& fields,
48 UnitCell<D>& unitCell) const
49 {
50 // Read header and check fields dimensions
51 int nMonomer;
52 bool isSymmetric;
53 readFieldHeader(in, nMonomer, unitCell, isSymmetric);
54 readMeshDimensions(in, mesh().dimensions());
55 checkAllocateFields(fields, nMonomer, mesh().dimensions());
56
57 // Allocate host arrays
59 allocateArrays(hostFields, nMonomer, mesh().size());
60
61 // Read data
62 Prdc::readRGridData(in, hostFields, nMonomer, mesh().dimensions());
63
64 // Copy device <- host
65 copyArrays(fields, hostFields);
66
67 // Return true iff the header contains a space group declaration
68 return isSymmetric;
69 }
70
71 /*
72 * Read the data section of an array of fields in r-grid format.
73 */
74 template <int D>
76 std::istream& in,
77 DArray< RField<D> >& fields,
78 int nMonomer) const
79 {
80 // Precondition: Check dimensions of fields
81 checkAllocateFields(fields, nMonomer, mesh().dimensions());
82
83 // Allocate host arrays
85 allocateArrays(hostFields, nMonomer, mesh().size());
86
87 // Read data section of file
88 Prdc::readRGridData(in, hostFields, nMonomer, mesh().dimensions());
89
90 // Copy device <- host
91 copyArrays(fields, hostFields);
92 }
93
94 /*
95 * Read a single field in r-grid format.
96 */
97 template <int D>
99 std::istream &in,
100 RField<D> & field,
101 UnitCell<D>& unitCell) const
102 {
103
104 // Read header and check field dimensions
105 int nMonomer;
106 bool isSymmetric;
107 readFieldHeader(in, nMonomer, unitCell, isSymmetric);
108 UTIL_CHECK(nMonomer == 1);
109 readMeshDimensions(in, mesh().dimensions());
110 checkAllocateField(field, mesh().dimensions());
111
112 // Allocate host field
113 HostDArray<cudaReal> hostField;
114 hostField.allocate(mesh().size());
115
116 // Read data section with one field
117 Prdc::readRGridData(in, hostField, mesh().dimensions());
118
119 // Copy device <- host
120 field = hostField;
121
122 // Return true iff the header contains a space group declaration
123 return isSymmetric;
124 }
125
126 /*
127 * Write an array of fields in r-grid format.
128 */
129 template <int D>
131 std::ostream &out,
132 DArray<RField<D> > const & fields,
133 UnitCell<D> const & unitCell,
134 bool writeHeader,
135 bool isSymmetric,
136 bool writeMeshSize) const
137 {
138 // Inspect fields array, check field dimensions
139 int nMonomer;
140 IntVec<D> meshDimensions;
141 inspectFields(fields, nMonomer, meshDimensions);
142 UTIL_CHECK(meshDimensions == mesh().dimensions());
143 int meshSize = mesh().size();
144 UTIL_CHECK(fields[0].capacity() == meshSize);
145
146 // Write header
147 if (writeHeader){
148 writeFieldHeader(out, nMonomer, unitCell, isSymmetric);
149 }
150 if (writeMeshSize){
151 writeMeshDimensions(out, meshDimensions);
152 }
153
154 // Copy field data to host container
155 DArray< HostDArray<cudaReal> > hostFields;
156 allocateArrays(hostFields, nMonomer, meshSize);
157 copyArrays(hostFields, fields);
158
159 // Write data section
160 Prdc::writeRGridData(out, hostFields, nMonomer, meshDimensions);
161 }
162
163 /*
164 * Write a single field in r-grid format.
165 */
166 template <int D>
168 std::ostream &out,
169 RField<D> const & field,
170 UnitCell<D> const & unitCell,
171 bool writeHeader,
172 bool isSymmetric) const
173 {
174 IntVec<D> meshDimensions = field.meshDimensions();
175 int meshSize = field.capacity();
176 UTIL_CHECK(meshDimensions == mesh().dimensions());
177 UTIL_CHECK(meshSize == mesh().size());
178
179 // Write header
180 if (writeHeader) {
181 writeFieldHeader(out, 1, unitCell, isSymmetric);
182 writeMeshDimensions(out, meshDimensions);
183 }
184
185 // Copy field (device) to hostField
186 HostDArray<cudaReal> hostField;
187 hostField.allocate(meshSize);
188 hostField = field;
189
190 // Write data from hostField
191 Prdc::writeRGridData(out, hostField, meshDimensions);
192 }
193
194 // Field IO in k-grid format
195
196 /*
197 * Read an array of fields in k-grid format
198 */
199 template <int D>
201 std::istream &in,
202 DArray<RFieldDft<D> >& fields,
203 UnitCell<D>& unitCell) const
204 {
205 // Read header and validate field mesh dimensions
206 int nMonomer;
207 bool isSymmetric;
208 readFieldHeader(in, nMonomer, unitCell, isSymmetric);
209 readMeshDimensions(in, mesh().dimensions());
210 checkAllocateFields(fields, nMonomer, mesh().dimensions());
211 IntVec<D> dftDimensions = fields[0].dftDimensions();
212 int capacity = fields[0].capacity();
213
214 // Allocate hostFields
216 allocateArrays(hostFields, nMonomer, capacity);
217
218 // Read data into hostFields
219 Prdc::readKGridData(in, hostFields, nMonomer, dftDimensions);
220
221 // Copy device <- host
222 copyArrays(fields, hostFields);
223 }
224
225 /*
226 * Write an array of fields in k-grid format
227 */
228 template <int D>
230 std::ostream &out,
231 DArray<RFieldDft<D> > const & fields,
232 UnitCell<D> const & unitCell,
233 bool isSymmetric) const
234 {
235 // Read header and validate field mesh dimensions
236 int nMonomer;
237 IntVec<D> meshDimensions;
238 inspectFields(fields, nMonomer, meshDimensions);
239 UTIL_CHECK(mesh().dimensions() == meshDimensions);
240 IntVec<D> dftDimensions = fields[0].dftDimensions();
241 int capacity = fields[0].capacity();
242
243 // Write header
244 writeFieldHeader(out, nMonomer, unitCell, isSymmetric);
245 writeMeshDimensions(out, meshDimensions);
246
247 // Copy data from device to hostFields
249 allocateArrays(hostFields, nMonomer, capacity);
250 copyArrays(hostFields, fields);
251
252 // Write data from hostFields
253 Prdc::writeKGridData(out, hostFields, nMonomer, dftDimensions);
254 }
255
256 /*
257 * Convert a single field from basis to k-grid format.
258 */
259 template <int D>
261 DArray<double> const & in,
262 RFieldDft<D>& out) const
263 {
265 UTIL_CHECK(out.isAllocated());
266 UTIL_CHECK(in.capacity() > 0);
267 UTIL_CHECK(out.meshDimensions() == mesh().dimensions());
268
269 // Allocate hostField
270 HostDArrayComplex hostField;
271 hostField.allocate(out.capacity());
272
273 // Convert basis to k-grid on hostField
274 Prdc::convertBasisToKGrid(in, hostField, basis(),
275 out.dftDimensions());
276
277 // Copy out (device) <- host
278 out = hostField;
279 }
280
281 /*
282 * Write an array of fields from k-grid to basis format.
283 */
284 template <int D>
286 RFieldDft<D> const & in,
287 DArray<double>& out,
288 bool checkSymmetry,
289 double epsilon) const
290 {
292 UTIL_CHECK(out.isAllocated());
293 UTIL_CHECK(in.meshDimensions() == mesh().dimensions());
294 UTIL_CHECK(out.capacity() > 0);
295
296 // Copy k-grid input to hostField
297 HostDArrayComplex hostField;
298 hostField.allocate(in.capacity());
299 hostField = in;
300
301 // Convert k-grid host field to basis format
302 Prdc::convertKGridToBasis(hostField, out, basis(),
303 in.dftDimensions(),
304 checkSymmetry, epsilon);
305 }
306
307 /*
308 * Test if an real field DFT has the declared space group symmetry.
309 */
310 template <int D>
312 RFieldDft<D> const & in,
313 double epsilon,
314 bool verbose) const
315 {
317 UTIL_CHECK(in.meshDimensions() == mesh().dimensions());
318
319 // Copy k-grid input to hostField
320 HostDArrayComplex hostField;
321 hostField.allocate(in.capacity());
322 hostField = in;
323
324 // Check symmetry of hostField
325 return Prdc::hasSymmetry(hostField, basis(), in.dftDimensions(),
326 epsilon, verbose);
327 }
328
329 /*
330 * Compare two fields in r-grid format, output report to Log file.
331 */
332 template <int D>
334 DArray< RField<D> > const & field2)
335 const
336 {
337 RFieldComparison<D> comparison;
338 comparison.compare(field1, field2);
339
340 Log::file() << "\n Real-space field comparison results"
341 << std::endl;
342 Log::file() << " Maximum Absolute Difference: "
343 << comparison.maxDiff() << std::endl;
344 Log::file() << " Root-Mean-Square Difference: "
345 << comparison.rmsDiff() << "\n" << std::endl;
346 }
347
348 /*
349 * Multiply a field in r-grid format by a constant factor.
350 */
351 template <int D>
353 RField<D> & field,
354 double factor) const
355 {
356 UTIL_CHECK(field.isAllocated());
357 VecOp::mulEqS(field, factor);
358 }
359
360 /*
361 * Replicate the unit cell for an array of r-grid fields.
362 */
363 template <int D>
365 std::ostream &out,
366 DArray< RField<D> > const & fields,
367 UnitCell<D> const & unitCell,
368 IntVec<D> const & replicas) const
369
370 {
371 // Inspect fields to obtain nMonomer and meshDimensions
372 int nMonomer;
373 IntVec<D> meshDimensions;
374 inspectFields(fields, nMonomer, meshDimensions);
375 UTIL_CHECK(meshDimensions == mesh().dimensions());
376 int capacity = fields[0].capacity();
377
378 // Copy k-grid input to hostField
379 DArray< HostDArray<cudaReal> > hostFields;
380 allocateArrays(hostFields, nMonomer, capacity);
381 copyArrays(hostFields, fields);
382
383 Prdc::replicateUnitCell(out, hostFields, meshDimensions,
384 unitCell, replicas);
385 }
386
387 /*
388 * Expand spatial dimension of an array of r-grid fields.
389 */
390 template <int D>
392 std::ostream &out,
393 DArray< RField<D> > const & fields,
394 UnitCell<D> const & unitCell,
395 int d,
396 DArray<int> const& newGridDimensions) const
397 {
398 // Inspect fields to obtain nMonomer and meshDimensions
399 int nMonomer;
400 IntVec<D> meshDimensions;
401 inspectFields(fields, nMonomer, meshDimensions);
402 UTIL_CHECK(meshDimensions == mesh().dimensions());
403 int capacity = fields[0].capacity();
404
405 // Copy k-grid input fields to hostFields
406 DArray< HostDArray<cudaReal> > hostFields;
407 allocateArrays(hostFields, nMonomer, capacity);
408 copyArrays(hostFields, fields);
409
410 Prdc::expandRGridDimension(out, hostFields, meshDimensions,
411 unitCell, d, newGridDimensions);
412 }
413
414}
415}
416#endif
double rmsDiff() const
Return the precomputed root-mean-squared difference.
double compare(FT const &a, FT const &b)
Compare individual fields.
double maxDiff() const
Return the precomputed maximum element-by-element difference.
Template for dynamic array stored in host CPU memory.
Definition HostDArray.h:43
An IntVec<D, T> is a D-component vector of elements of integer type T.
Definition IntVec.h:27
bool isAllocated() const
Return true if the FftwDArray has been allocated, false otherwise.
Definition FftwDArray.h:102
Comparator for fields in real-space (r-grid) format.
Fourier transform of a real field on an FFT mesh.
IntVec< D > const & dftDimensions() const
Return vector of dft (Fourier) grid dimensions by constant reference.
IntVec< D > const & meshDimensions() const
Return vector of spatial mesh dimensions by constant reference.
Field of real double precision values on an FFT mesh.
Definition cpu/RField.h:29
const IntVec< D > & meshDimensions() const
Return mesh dimensions by constant reference.
Definition cpu/RField.h:112
Base template for UnitCell<D> classes, D=1, 2 or 3.
Definition UnitCell.h:56
bool readFieldRGrid(std::istream &in, RField< D > &field, UnitCell< D > &unitCell) const override
Read a single RField (field on an r-space grid) from a stream.
void convertKGridToBasis(RFieldDft< D > const &in, DArray< double > &out, bool checkSymmetry=true, double epsilon=1.0e-8) const override
Convert a field from Fourier (k-grid) to symmetrized basis form.
void convertBasisToKGrid(DArray< double > const &components, RFieldDft< D > &dft) const override
Convert a field from symmetrized basis to Fourier grid (k-grid).
Basis< D > const & basis() const
Get the associated Basis by const reference.
void compareFieldsRGrid(DArray< RField< D > > const &field1, DArray< RField< D > > const &field2) const override
Compare two fields in r-grid format, output a report.
bool hasSymmetry(RFieldDft< D > const &in, double epsilon=1.0e-8, bool verbose=true) const override
Check if a k-grid field has the declared space group symmetry.
void writeFieldsRGrid(std::ostream &out, DArray< RField< D > > const &fields, UnitCell< D > const &unitCell, bool writeHeader=true, bool isSymmetric=true, bool writeMeshSize=true) const override
Write array of RField objects (fields on r-space grid) to a stream.
void expandRGridDimension(std::ostream &out, DArray< RField< D > > const &fields, UnitCell< D > const &unitCell, int d, DArray< int > const &newGridDimensions) const override
Expand spatial dimension of an array of r-grid fields.
Mesh< D > const & mesh() const
Get spatial discretization mesh by const reference.
void writeFieldsKGrid(std::ostream &out, DArray< RFieldDft< D > > const &fields, UnitCell< D > const &unitCell, bool isSymmetric=true) const override
Write array of RFieldDft objects (k-space fields) to file.
void readFieldsKGrid(std::istream &in, DArray< RFieldDft< D > > &fields, UnitCell< D > &unitCell) const override
Read array of RFieldDft objects (k-space fields) from a stream.
bool readFieldsRGrid(std::istream &in, DArray< RField< D > > &fields, UnitCell< D > &unitCell) const override
Read array of RField objects (r-grid fields) from a stream.
void readFieldsRGridData(std::istream &in, DArray< RField< D > > &fields, int nMonomer) const override
Read data for an array of r-grid fields, with no header section.
void replicateUnitCell(std::ostream &out, DArray< RField< D > > const &fields, UnitCell< D > const &unitCell, IntVec< D > const &replicas) const override
Write r-grid fields in a replicated unit cell to std::ostream.
void writeFieldRGrid(std::ostream &out, RField< D > const &field, UnitCell< D > const &unitCell, bool writeHeader=true, bool isSymmetric=true) const override
Write a single RField (field on an r-space grid) to a stream.
void scaleFieldRGrid(RField< D > &field, double factor) const override
Rescale a single r-grid field by a scalar factor.
HostDArray containing cudaComplex elements.
int capacity() const
Return allocated size.
Definition Array.h:159
Dynamically allocatable contiguous array template.
Definition DArray.h:32
void allocate(int capacity)
Allocate the underlying C array.
Definition DArray.h:199
bool isAllocated() const
Return true if this DArray has been allocated, false otherwise.
Definition DArray.h:247
static std::ostream & file()
Get log ostream by reference.
Definition Log.cpp:59
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
void replicateUnitCell(IntVec< 1 > const &replicas, UnitCell< 1 > const &cellIn, UnitCell< 1 > &cellOut)
Create a replicated UnitCell<1>.
void writeFieldHeader(std::ostream &out, int ver1, int ver2, UnitCell< D > const &cell, std::string const &groupName, int nMonomer)
Write common part of field header (fortran PSCF format).
void readFieldHeader(std::istream &in, int &ver1, int &ver2, UnitCell< D > &cell, std::string &groupName, int &nMonomer)
Read common part of field header (fortran PSCF format).
void writeKGridData(std::ostream &out, DArray< AT > const &fields, int nMonomer, IntVec< D > const &dftDimensions)
Write data for array of k-grid fields, with no header section.
bool hasSymmetry(AT const &in, Basis< D > const &basis, IntVec< D > const &dftDimensions, double epsilon=1.0e-8, bool verbose=true)
Check if a k-grid field has the declared space group symmetry.
void inspectFields(DArray< FT > const &fields, int &nMonomer, IntVec< D > &dimensions)
Inspect dimensions of a DArray of fields, each of type FT.
void readKGridData(std::istream &in, DArray< AT > &fields, int nMonomer, IntVec< D > const &dftDimensions)
Read data for array of k-grid fields, with no header section.
void convertBasisToKGrid(DArray< double > const &components, AT &dft, Basis< D > const &basis, IntVec< D > const &dftDimensions)
Convert a real field from symmetrized basis to Fourier grid.
void writeMeshDimensions(std::ostream &out, IntVec< D > const &meshDimensions)
Write mesh dimensions to a field file header.
void expandRGridDimension(std::ostream &out, DArray< AT > const &fields, IntVec< D > const &meshDimensions, UnitCell< D > const &unitCell, int d, DArray< int > newGridDimensions)
Expand the dimensionality of space from D to d.
void writeRGridData(std::ostream &out, DArray< AT > const &fields, int nMonomer, IntVec< D > const &dimensions)
Write data for array of r-grid fields, with no header section.
void convertKGridToBasis(AT const &in, DArray< double > &out, Basis< D > const &basis, IntVec< D > const &dftDimensions, bool checkSymmetry=true, double epsilon=1.0e-8)
Convert a real field from Fourier grid to symmetrized basis.
void readRGridData(std::istream &in, DArray< AT > &fields, int nMonomer, IntVec< D > const &dimensions)
Read data for array of r-grid fields, with no header section.
void readMeshDimensions(std::istream &in, IntVec< D > const &meshDimensions)
Read mesh dimensions from a field file header.
void checkAllocateFields(DArray< FT > &fields, int nMonomer, IntVec< D > const &dimensions)
Check allocation of an array of fields, allocate if necessary.
void checkAllocateField(FT &field, IntVec< D > const &dimensions)
Check allocation of a single field, allocate if necessary.
void mulEqS(DeviceArray< cudaReal > &a, const cudaReal b, const int beginIdA, const int n)
Vector multiplication in-place, a[i] *= b, kernel wrapper (cudaReal).
Definition VecOp.cu:1918
Fields, FFTs, and utilities for periodic boundary conditions (CUDA)
Definition Reduce.cpp:14
Periodic fields and crystallography.
Definition CField.cpp:11
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.
Definition param_pc.dox:1