1 /* Copyright (C) 1999-2020 Massachusetts Institute of Technology.
2 *
3 * This program is free software; you can redistribute it and/or modify
4 * it under the terms of the GNU General Public License as published by
5 * the Free Software Foundation; either version 2 of the License, or
6 * (at your option) any later version.
7 *
8 * This program is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * GNU General Public License for more details.
12 *
13 * You should have received a copy of the GNU General Public License
14 * along with this program; if not, write to the Free Software
15 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16 */
17
18 #include <stdlib.h>
19 #include <stdio.h>
20 #include <string.h>
21 #include <math.h>
22
23 #include "config.h"
24
25 #include <check.h>
26 #include <matrices.h>
27
28 #include "matrixio.h"
29
30 #define TWOPI 6.2831853071795864769252867665590057683943388
31 #define MAX2(a,b) ((a) > (b) ? (a) : (b))
32
33 /* note that kvector here is given in the reciprocal basis
34 ...data_id should be of length at 2*num_components */
fieldio_write_complex_field(scalar_complex * field,int rank,const int dims[3],const int local_dims[3],const int start[3],int which_component,int num_components,const real * kvector,matrixio_id file_id,int append,matrixio_id data_id[])35 void fieldio_write_complex_field(scalar_complex *field,
36 int rank,
37 const int dims[3],
38 const int local_dims[3],
39 const int start[3],
40 int which_component, int num_components,
41 const real *kvector,
42 matrixio_id file_id,
43 int append,
44 matrixio_id data_id[])
45 {
46 int i, j, k, component, ri_part;
47
48 rank = dims[2] == 1 ? (dims[1] == 1 ? 1 : 2) : 3;
49
50 if (kvector) {
51 real s[3]; /* the step size between grid points dotted with k */
52 scalar_complex *phasex, *phasey, *phasez;
53
54 for (i = 0; i < 3; ++i)
55 s[i] = TWOPI * kvector[i] / dims[i];
56
57 /* cache exp(ikx) along each of the directions, for speed */
58 CHK_MALLOC(phasex, scalar_complex, local_dims[0]);
59 CHK_MALLOC(phasey, scalar_complex, local_dims[1]);
60 CHK_MALLOC(phasez, scalar_complex, local_dims[2]);
61 for (i = 0; i < local_dims[0]; ++i) {
62 real phase = s[0] * (i + start[0]);
63 phasex[i].re = cos(phase);
64 phasex[i].im = sin(phase);
65 }
66 for (j = 0; j < local_dims[1]; ++j) {
67 real phase = s[1] * (j + start[1]);
68 phasey[j].re = cos(phase);
69 phasey[j].im = sin(phase);
70 }
71 for (k = 0; k < local_dims[2]; ++k) {
72 real phase = s[2] * (k + start[2]);
73 phasez[k].re = cos(phase);
74 phasez[k].im = sin(phase);
75 }
76
77 /* Now, multiply field by exp(i k*r): */
78 for (i = 0; i < local_dims[0]; ++i) {
79 scalar_complex px = phasex[i];
80
81 for (j = 0; j < local_dims[1]; ++j) {
82 scalar_complex py;
83 real re = phasey[j].re, im = phasey[j].im;
84 py.re = px.re * re - px.im * im;
85 py.im = px.re * im + px.im * re;
86
87 for (k = 0; k < local_dims[2]; ++k) {
88 int ijk = ((i*local_dims[1] + j)*local_dims[2] + k)*3;
89 real p_re, p_im;
90 real re = phasez[k].re, im = phasez[k].im;
91
92 p_re = py.re * re - py.im * im;
93 p_im = py.re * im + py.im * re;
94
95 for (component = 0; component < 3; ++component) {
96 int ijkc = ijk + component;
97 re = field[ijkc].re; im = field[ijkc].im;
98 field[ijkc].re = re * p_re - im * p_im;
99 field[ijkc].im = im * p_re + re * p_im;
100 }
101 }
102 }
103 }
104
105 free(phasez);
106 free(phasey);
107 free(phasex);
108 }
109
110 /* write hyperslabs for each field component: */
111 for (component = 0; component < num_components; ++component)
112 if (component == which_component ||
113 which_component < 0)
114 for (ri_part = 0; ri_part < 2; ++ri_part) {
115 char name[] = "x.i";
116 name[0] = (num_components == 1 ? 'c' : 'x') + component;
117 name[2] = ri_part ? 'i' : 'r';
118
119 if (!append)
120 data_id[component*2 + ri_part] =
121 matrixio_create_dataset(file_id, name, NULL,
122 rank, dims);
123
124 matrixio_write_real_data(
125 data_id[component*2 + ri_part], local_dims, start,
126 2 * num_components,
127 ri_part ? &field[component].im
128 : &field[component].re);
129 }
130 }
131
fieldio_write_real_vals(real * vals,int rank,const int dims[3],const int local_dims[3],const int start[3],matrixio_id file_id,int append,const char * dataname,matrixio_id * data_id)132 void fieldio_write_real_vals(real *vals,
133 int rank,
134 const int dims[3],
135 const int local_dims[3],
136 const int start[3],
137 matrixio_id file_id,
138 int append,
139 const char *dataname,
140 matrixio_id *data_id)
141 {
142 rank = dims[2] == 1 ? (dims[1] == 1 ? 1 : 2) : 3;
143
144 if (!append || data_id->id < 0)
145 *data_id = matrixio_create_dataset(file_id, dataname,
146 NULL, rank,dims);
147
148 matrixio_write_real_data(*data_id,local_dims,start,1,vals);
149 }
150