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