oILAB
Loading...
Searching...
No Matches
Diff.h
Go to the documentation of this file.
1//
2// Created by Nikhil Chandra Admal on 7/14/23.
3//
4
5#ifndef OILAB_DIFF_H
6#define OILAB_DIFF_H
7
8#include <Operator.h>
9#include <LatticeModule.h>
10#include <unsupported/Eigen/CXX11/Tensor>
11#include <FFT.h>
12
13
14namespace gbLAB {
15 template<int dim>
16 class Diff : public Operator<Diff<dim>,dim>
17 {
18 public:
19 using dcomplex= std::complex<double>;
20 using Operator<Diff<dim>,dim>::n, Operator<Diff<dim>,dim>::L;
21 const Eigen::array<Eigen::Index,dim> d;
22
23 explicit Diff(const Eigen::array<Eigen::Index,dim>& d_,
24 const Eigen::Matrix<double,dim,dim>& A,
25 const Eigen::array<Eigen::Index,dim>& n_) :
26 Operator<Diff<dim>,dim>(A,n_),
27 d(d_)
28 {
29 for(const auto& order : d)
30 assert(order >= 0);
31
32 }
33
34 template<int dm=dim>
35 typename std::enable_if<dm==1,void>::type
36 perform_op(const double* x_in, double* y_out) const
37 {
38 Eigen::TensorMap<const Eigen::Tensor<double,dm>> xReal(x_in,n);
39 const Eigen::Tensor<dcomplex,dm> x=xReal.template cast<dcomplex>();
40 Eigen::TensorMap<Eigen::Tensor<double,dm>> y(y_out,n);
41
42 // if d=0 y_out= x_in and return
43 int totalOrder= 0;
44 for(const auto& order : d)
45 {
46 totalOrder= totalOrder + order;
47 }
48 if (totalOrder == 0) {
49 y= xReal;
50 return;
51 }
52
53 // Compute y = Lx using FFT
54 Eigen::Tensor<dcomplex,dm> xhat(n);
55 xhat.setZero();
56 FFT::fft(x,xhat);
57
58 Eigen::Tensor<dcomplex,dm> d2fhat(n);
59 d2fhat.setZero();
60
61 // only part which is dimension dependent
62 for (int i = 0; i < n[0]; ++i) {
63 ReciprocalLatticeVector<dm> r(L);
64 dcomplex factor(1,0);
65 for (int k= 0; k<dim; ++k) {
66 if (d[k] == 0) continue;
67 else if (d[k] % 2 == 0)
68 r << (i <= n[0] / 2 ? i : i - n[0]);
69 else
70 r << (i == n[0] / 2 ? 0 : -n[0] * (i / (n[0] / 2)) + i);
71 factor = factor * std::pow(-2.0 * M_PI *
72 dcomplex(0, 1) *
73 r.cartesian()(k),
74 d[k]);
75 }
76 d2fhat(i)= xhat(i) * factor;
77 }
78
79 // Laplacian Lf
80 Eigen::Tensor<dcomplex,dm> Lf(n);
81 Lf.setZero();
82 FFT::ifft(d2fhat,Lf);
83 y= Lf.real();
84 }
85
86 template<int dm=dim>
87 typename std::enable_if<dm==2,void>::type
88 perform_op(const double* x_in, double* y_out) const
89 {
90 Eigen::TensorMap<const Eigen::Tensor<double,2>> xReal(x_in,n);
91 const Eigen::Tensor<dcomplex,2> x=xReal.cast<dcomplex>();
92 Eigen::TensorMap<Eigen::Tensor<double,2>> y(y_out,n);
93
94 // if d=0 y_out= x_in and return
95 int totalOrder= 0;
96 for(const auto& order : d)
97 {
98 totalOrder= totalOrder + order;
99 }
100 if (totalOrder == 0) {
101 y= xReal;
102 return;
103 }
104
105 // Compute y = Lx using FFT
106 Eigen::Tensor<dcomplex,2> xhat(n);
107 xhat.setZero();
108 FFT::fft(x,xhat);
109
110 Eigen::Tensor<dcomplex,2> d2fhat(n);
111 d2fhat.setZero();
112
113 // only part which is dimension dependent
114 for (int i = 0; i < n[0]; ++i) {
115 for (int j = 0; j < n[1]; ++j) {
116 ReciprocalLatticeVector<2> r(L);
117 dcomplex factor(1,0);
118 for (int k= 0; k<dim; ++k) {
119 if (d[k] == 0) continue;
120 else if (d[k] % 2 == 0)
121 r << (i <= n[0] / 2 ? i : i - n[0]),
122 (j <= n[1] / 2 ? j : j - n[1]);
123 else
124 r << (i == n[0] / 2 ? 0 : -n[0] * (i / (n[0] / 2)) + i),
125 (j == n[1] / 2 ? 0 : -n[1] * (j / (n[1] / 2)) + j);
126 factor = factor * std::pow(-2.0 * M_PI *
127 dcomplex(0, 1) *
128 r.cartesian()(k),
129 d[k]);
130 }
131 d2fhat(i,j)= xhat(i,j) * factor;
132 }
133 }
134
135 // Laplacian Lf
136 Eigen::Tensor<dcomplex,2> Lf(n);
137 Lf.setZero();
138 FFT::ifft(d2fhat,Lf);
139 y= Lf.real();
140 }
141
142 template<int dm=dim>
143 typename std::enable_if<dm==3,void>::type
144 perform_op(const double* x_in, double* y_out) const
145 {
146 Eigen::TensorMap<const Eigen::Tensor<double,dm>> xReal(x_in,n);
147 const Eigen::Tensor<dcomplex,dm> x=xReal.template cast<dcomplex>();
148 Eigen::TensorMap<Eigen::Tensor<double,dm>> y(y_out,n);
149
150 // if d=0 y_out= x_in and return
151 int totalOrder= 0;
152 for(const auto& order : d)
153 {
154 totalOrder= totalOrder + order;
155 }
156 if (totalOrder == 0) {
157 y= xReal;
158 return;
159 }
160
161 // Compute y = Lx using FFT
162 Eigen::Tensor<dcomplex,dm> xhat(n);
163 xhat.setZero();
164 FFT::fft(x,xhat);
165
166 Eigen::Tensor<dcomplex,dm> d2fhat(n);
167 d2fhat.setZero();
168
169 // only part which is dimension dependent
170 for (int i = 0; i < n[0]; ++i) {
171 for (int j = 0; j < n[1]; ++j) {
172 for (int k = 0; k < n[2]; ++k) {
173 ReciprocalLatticeVector<dm> r(L);
174 dcomplex factor(1, 0);
175 for (int l = 0; l < dm; ++l) {
176 if (d[l] == 0) continue;
177 else if (d[l] % 2 == 0)
178 r << (i <= n[0] / 2 ? i : i - n[0]),
179 (j <= n[1] / 2 ? j : j - n[1]),
180 (k <= n[2] / 2 ? k : k - n[2]);
181 else
182 r << (i == n[0] / 2 ? 0 : -n[0] * (i / (n[0] / 2)) + i),
183 (j == n[1] / 2 ? 0 : -n[1] * (j / (n[1] / 2)) + j),
184 (k == n[2] / 2 ? 0 : -n[2] * (k / (n[2] / 2)) + k);
185 factor = factor * std::pow(-2.0 * M_PI *
186 dcomplex(0, 1) *
187 r.cartesian()(l),
188 d[l]);
189 }
190 d2fhat(i,j,k) = xhat(i,j,k) * factor;
191 }
192 }
193 }
194
195 // Laplacian Lf
196 Eigen::Tensor<dcomplex,dm> Lf(n);
197 Lf.setZero();
198 FFT::ifft(d2fhat,Lf);
199 y= Lf.real();
200 }
201 };
202}
203
204#endif //OILAB_DIFF_H
FFT::ifft
static void ifft(const Eigen::Tensor< dcomplex, 3 > &in, Eigen::Tensor< dcomplex, 3 > &out)
Definition FFT.h:54
FFT::fft
static void fft(const Eigen::Tensor< dcomplex, 3 > &in, Eigen::Tensor< dcomplex, 3 > &out)
Definition FFT.h:17
gbLAB::Diff::perform_op
std::enable_if< dm==1, void >::type perform_op(const double *x_in, double *y_out) const
Definition Diff.h:36
gbLAB::Diff::perform_op
std::enable_if< dm==2, void >::type perform_op(const double *x_in, double *y_out) const
Definition Diff.h:88
gbLAB::Diff::d
const Eigen::array< Eigen::Index, dim > d
Definition Diff.h:21
gbLAB::Diff::dcomplex
std::complex< double > dcomplex
Definition Diff.h:19
gbLAB::Diff::Diff
Diff(const Eigen::array< Eigen::Index, dim > &d_, const Eigen::Matrix< double, dim, dim > &A, const Eigen::array< Eigen::Index, dim > &n_)
Definition Diff.h:23
gbLAB::Diff::perform_op
std::enable_if< dm==3, void >::type perform_op(const double *x_in, double *y_out) const
Definition Diff.h:144
gbLAB::Operator< Diff< dim >, dim >::n
const Eigen::array< Eigen::Index, dim > n
Definition Operator.h:18