MADNESS  0.10.1
DKops.h
Go to the documentation of this file.
1 #ifndef MADNESS_APPS_MOLDFT_DKOPS_H_INCLUDED
2 #define MADNESS_APPS_MOLDFT_DKOPS_H_INCLUDED
3 
4 
5 #include <madness/mra/mra.h>
6 #include <madness/constants.h>
7 #include <madness/mra/operator.h>
8 #include <iostream>
9 #include <fstream>
10 #include <exception>
11 
12 //#include "../../../../mpfrc++-3/mpreal.h"
13 
14 using namespace madness;
15 
16 static const double opthresh = 1e-16;
17 
18 double q(double t){
19  return exp(-t);
20 }
21 
22 double w(double t, double eps){
23  double c = 137.0359895;
24  double PI = constants::pi;
25  return 1.0/(c*sqrt(PI))*exp(-1.0/2.0*t-c*c*exp(-t))-(1.0+eps/(c*c))*exp((2.0*eps+(eps*eps)/(c*c))*exp(-t)-t)*(erfc((c+eps/c)*(exp(-t/2.0)))-2.0);
26 }
27 
28 double Ebark(double t, double eps){
29  double R = 1e-8;
30  return pow(1.0/(2.0*q(t)), 1.5)*w(t, eps)*exp(-1.0/(4.0*q(t))*R*R);
31 }
32 
33 real_convolution_3d Ebar(World& world, double eps){
34  Vector<double,3> args = vec(0.0,0.0,0.0);
35 
36 
37  std::vector<double> cvec;
38  std::vector<double> tvec;
39 
40  double fac = pow(2.0*constants::pi,-1.5);
41  double dt = 1.0/8.0;
42  double tval = -10.0;
43  while(tval < 100.0){
44 
45  if(Ebark(tval, eps) > opthresh){
46  //if(world.rank()==0) print(" Ebark(",tval,",",eps," = ",Ebark(tval,eps));
47  cvec.push_back(dt/pow(2.0*q(tval),1.5)*w(tval,eps)*fac);
48  tvec.push_back(1.0/(4.0*q(tval)));
49  }
50  //else{
51  //if(world.rank()==0) print("Ebark(",tval,",",eps," = ",Ebark(tval,eps));
52 
53  //}
54  tval = tval + dt;
55  }
56 
57  Tensor<double> ctensor(cvec.size());
58  Tensor<double> ttensor(tvec.size());
59  for(int i = 0; i < cvec.size(); i++){
60  ctensor[i] = cvec[i];
61  ttensor[i] = tvec[i];
62  }
63 
64  int n = cvec.size();
65  //if(world.rank()==0) print("Made an Ebar! n = ",n);
66  //if(world.rank()==0) print("Made an Ebar! Here's what's inside:\n\nc:\n",ctensor,"\nt:\n",ttensor);
67 
68  //test ebar?
69  /*
70  if(world.rank()==0) print("Testing Ebar. Value should be ~22.680321967196");
71  double testvalue = 0.0;
72  for(int i = 0; i < cvec.size(); i++){
73  testvalue += ctensor[i]*exp(-1e-2*ttensor[i]);
74  }
75  if(world.rank()==0) print("And the result is: ", testvalue);
76  */
77 
78  //return real_convolution_3d(world, args, ctensor, ttensor);
79  return real_convolution_3d(world, ctensor, ttensor);
80 }
81 
82 //real_convolution_3d Ebar_fixed(World& world){
83 // //Tensor<double> c(305l), t(305l);
84 // std::vector<double> c(365), t(365);
85 // Vector<double,3> args = vec(0.0,0.0,0.0);
86 // #include "RelCoeffs/Ebar_coeffs.dat"
87 // int n = 305;
88 //
89 // for(int i=0; i < n; i++){
90 // if(c[i]*exp(-t[i]*1e-16) < opthresh){
91 // c.erase(c.begin()+i);
92 // t.erase(t.begin()+i);
93 // i--;
94 // n--;
95 // }
96 // }
97 //
98 // //print("n = ", n);
99 //
100 // Tensor<double> ctens(n), ttens(n);
101 // for(int i = 0; i < n; i++){
102 // ctens[i] = c[i];
103 // ttens[i] = t[i];
104 // }
105 // //if(world.rank()==0) print("Made a Ebar (fixed)! n = ", n);
106 // //if(world.rank()==0) print("Made a Pbar! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
107 // //return real_convolution_3d(world, args, ctens, ttens);
108 // return real_convolution_3d(world, ctens, ttens);
109 //}
110 
111 real_convolution_3d Pbar(World& world){
112  //Tensor<double> c(305l), t(305l);
113  std::vector<double> c;
114  std::vector<double> t;
115  Vector<double,3> args = vec(0.0,0.0,0.0);
116 
117  std::ifstream inf("/gpfs/home/jscanderson/DKproject/Pbar_t.csv");
118  if(!inf){
119  if(world.rank() == 0) std::cerr << "Unable to open Pbar_t.csv" << std::endl;
120  exit(1);
121  }
122  std::string strInput;
123  getline(inf, strInput);
124  while(inf){
125  t.push_back(std::stod(strInput));
126  getline(inf, strInput);
127  }
128  inf.close();
129  //if(world.rank()==0) print("Done Reading Pbar_t");
130 
131  inf.open("/gpfs/home/jscanderson/DKproject/Pbar_c.csv");
132  if(!inf){
133  if(world.rank() == 0) std::cerr << "Unable to open Pbar_c.csv" << std::endl;
134  exit(1);
135  }
136  getline(inf, strInput);
137  while(inf){
138  //if(world.rank()==0) print("a", strInput, "b");
139  c.push_back(std::stod(strInput));
140  getline(inf, strInput);
141  }
142  inf.close();
143  //if(world.rank()==0) print("Done Reading Pbar_c");
144 
145 
146  int n = c.size();
147  /*
148  for(int i=0; i < n; i++){
149  if(c[i]*exp(-t[i]*1e-30) < opthresh){
150  c.erase(c.begin()+i);
151  t.erase(t.begin()+i);
152  i--;
153  n--;
154  }
155  }
156  */
157  //print("n = ", n);
158 
159  Tensor<double> ctens(n), ttens(n);
160  for(int i = 0; i < n; i++){
161  ctens[i] = c[i];
162  ttens[i] = t[i];
163  }
164  //if(world.rank()==0) print("Made a Pbar! n = ", n);
165  //if(world.rank()==0) print("Made a Pbar! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
166  //return real_convolution_3d(world, args, ctens, ttens);
167  return real_convolution_3d(world, ctens, ttens);
168 }
169 
170 real_convolution_3d Tbar(World& world){
171  //Tensor<double> c(305l), t(305l);
172  std::vector<double> c;
173  std::vector<double> t;
174  Vector<double,3> args = vec(0.0,0.0,0.0);
175 
176  std::ifstream inf("/gpfs/home/jscanderson/DKproject/Tbar_t.csv");
177  if(!inf){
178  if(world.rank() == 0) std::cerr << "Unable to open Tbar_t.csv" << std::endl;
179  exit(1);
180  }
181  std::string strInput;
182  getline(inf, strInput);
183  while(inf){
184  t.push_back(std::stod(strInput));
185  getline(inf, strInput);
186  }
187  inf.close();
188  //if(world.rank()==0) print("Done Reading Tbar_t");
189 
190  inf.open("/gpfs/home/jscanderson/DKproject/Tbar_c.csv");
191  if(!inf){
192  if(world.rank() == 0) std::cerr << "Unable to open Tbar_c.csv" << std::endl;
193  exit(1);
194  }
195  getline(inf, strInput);
196  while(inf){
197  //if(world.rank()==0) print("a", strInput, "b");
198  c.push_back(std::stod(strInput));
199  getline(inf, strInput);
200  }
201  inf.close();
202  //if(world.rank()==0) print("Done Reading Tbar_c");
203 
204 
205  int n = c.size();
206  /*
207  for(int i=0; i < n; i++){
208  if(c[i]*exp(-t[i]*1e-30) < opthresh){
209  c.erase(c.begin()+i);
210  t.erase(t.begin()+i);
211  i--;
212  n--;
213  }
214  }
215  */
216  //print("n = ", n);
217 
218  Tensor<double> ctens(n), ttens(n);
219  for(int i = 0; i < n; i++){
220  ctens[i] = c[i];
221  ttens[i] = t[i];
222  }
223  //if(world.rank()==0) print("Made a Tbar! n = ", n);
224  //if(world.rank()==0) print("Made a Pbar! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
225  //return real_convolution_3d(world, args, ctens, ttens);
226  return real_convolution_3d(world, ctens, ttens);
227 }
228 
229 
230 real_convolution_3d A(World& world){
231  //Tensor<double> c(301l), t(301l);
232  std::vector<double> c, t;
233  Vector<double,3> args = vec(0.0,0.0,0.0);
234 
235  std::ifstream inf("/gpfs/home/jscanderson/DKproject/A_t.csv");
236  if(!inf){
237  if(world.rank() == 0) std::cerr << "Unable to open A_t.csv" << std::endl;
238  exit(1);
239  }
240  std::string strInput;
241  getline(inf, strInput);
242  while(inf){
243  t.push_back(std::stod(strInput));
244  getline(inf, strInput);
245  }
246  inf.close();
247  //if(world.rank()==0) print("Done Reading A_t");
248 
249  inf.open("/gpfs/home/jscanderson/DKproject/A_c.csv");
250  if(!inf){
251  if(world.rank() == 0) std::cerr << "Unable to open A_c.csv" << std::endl;
252  exit(1);
253  }
254  getline(inf, strInput);
255  while(inf){
256  c.push_back(std::stod(strInput));
257  getline(inf, strInput);
258  }
259  inf.close();
260  //if(world.rank()==0) print("Done Reading A_c");
261 
262  int n = c.size();
263  /*
264  for(int i=0; i < n; i++){
265  if(c[i]*exp(-t[i]*1e-30) < opthresh){
266  c.erase(c.begin()+i);
267  t.erase(t.begin()+i);
268  i--;
269  n--;
270  }
271  }
272  */
273  double fac = pow(2.0*constants::pi,1.5);
274  Tensor<double> ctens(n), ttens(n);
275  for(int i = 0; i < n; i++){
276  ctens[i] = c[i]*fac;
277  ttens[i] = t[i];
278  }
279  //if(world.rank()==0) print("Made an A!, n = ", n);
280  //if(world.rank()==0) print("Made an A! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
281  //return real_convolution_3d(world, args, ctens, ttens);
282  return real_convolution_3d(world, ctens, ttens);
283 }
284 
285 //Same as A(World& world) above but adds 1/sqrt(2) to the largest coefficient
286 real_convolution_3d A2(World& world){
287  //Tensor<double> c(301l), t(301l);
288  std::vector<double> c, t;
289  Vector<double,3> args = vec(0.0,0.0,0.0);
290 
291  std::ifstream inf("/gpfs/home/jscanderson/DKproject/A_t.csv");
292  if(!inf){
293  if(world.rank() == 0) std::cerr << "Unable to open A_t.csv" << std::endl;
294  exit(1);
295  }
296  std::string strInput;
297  getline(inf, strInput);
298  while(inf){
299  t.push_back(std::stod(strInput));
300  getline(inf, strInput);
301  }
302  inf.close();
303  //if(world.rank()==0) print("Done Reading A_t");
304 
305  inf.open("/gpfs/home/jscanderson/DKproject/A_c.csv");
306  if(!inf){
307  if(world.rank() == 0) std::cerr << "Unable to open A_c.csv" << std::endl;
308  exit(1);
309  }
310  getline(inf, strInput);
311  while(inf){
312  c.push_back(std::stod(strInput));
313  getline(inf, strInput);
314  }
315  inf.close();
316  //if(world.rank()==0) print("Done Reading A_c");
317 
318  int n = c.size();
319  /*
320  for(int i=0; i < n; i++){
321  if(c[i]*exp(-t[i]*1e-30) < opthresh){
322  c.erase(c.begin()+i);
323  t.erase(t.begin()+i);
324  i--;
325  n--;
326  }
327  }
328  */
329  double fac = pow(2.0*constants::pi,1.5);
330  Tensor<double> ctens(n), ttens(n);
331  for(int i = 0; i < n; i++){
332  ctens[i] = c[i]*fac;
333  ttens[i] = t[i];
334  }
335  c[n-1] += 1/sqrt(2)*fac*std::pow(t[n-1]*constants::pi,1.5); //largest exponent is at the end, so add 1/sqrt(2) to the coefficient there
336  //if(world.rank()==0) print("Made an A!, n = ", n);
337  //if(world.rank()==0) print("Made an A! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
338  //return real_convolution_3d(world, args, ctens, ttens);
339  return real_convolution_3d(world, ctens, ttens);
340 }
341 
342 real_convolution_3d PbarA(World& world){
343  //Tensor<double> c(441l), t(441l);
344  std::vector<double> c, t;
345  Vector<double,3> args = vec(0.0,0.0,0.0);
346 
347  std::ifstream inf("/gpfs/home/jscanderson/DKproject/PbarA_t.csv");
348  if(!inf){
349  if(world.rank() == 0) std::cerr << "Unable to open PbarA_t.csv" << std::endl;
350  exit(1);
351  }
352  std::string strInput;
353  getline(inf, strInput);
354  while(inf){
355  t.push_back(std::stod(strInput));
356  getline(inf, strInput);
357  }
358  inf.close();
359  //if(world.rank()==0) print("Done Reading PbarA_t");
360 
361  inf.open("/gpfs/home/jscanderson/DKproject/PbarA_c.csv");
362  if(!inf){
363  if(world.rank() == 0) std::cerr << "Unable to open PbarA_c.csv" << std::endl;
364  exit(1);
365  }
366  getline(inf, strInput);
367  while(inf){
368  c.push_back(std::stod(strInput));
369  getline(inf, strInput);
370  }
371  inf.close();
372  //if(world.rank()==0) print("Done Reading PbarA_c");
373 
374  int n = c.size();
375  /*
376  for(int i=0; i < n; i++){
377  if(c[i]*exp(-t[i]*1e-30) < opthresh){
378  c.erase(c.begin()+i);
379  t.erase(t.begin()+i);
380  i--;
381  n--;
382  }
383  }
384  */
385  double fac = pow(2.0*constants::pi,1.5);
386  Tensor<double> ctens(n), ttens(n);
387  for(int i = 0; i < n; i++){
388  ctens[i] = c[i]*fac;
389  ttens[i] = t[i];
390  }
391  //if(world.rank()==0) print("Made a PbarA!, n = ", n);
392  //if(world.rank()==0) print("Made a PbarA! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
393  //return real_convolution_3d(world, args, ctens, ttens);
394  return real_convolution_3d(world, ctens, ttens);
395 }
396 
397 std::vector<real_convolution_3d_ptr> gradPbarA(World& world){
398  const double pi = constants::pi;
399  std::vector<double> c, t;
400 
401  std::ifstream inf("/gpfs/home/jscanderson/DKproject/PbarA_t.csv");
402  if(!inf){
403  if(world.rank() == 0) std::cerr << "Unable to open PbarA_t.csv" << std::endl;
404  exit(1);
405  }
406  std::string strInput;
407  getline(inf, strInput);
408  while(inf){
409  t.push_back(std::stod(strInput));
410  getline(inf, strInput);
411  }
412  inf.close();
413  //if(world.rank()==0) print("Done Reading PbarA_t");
414 
415  inf.open("/gpfs/home/jscanderson/DKproject/PbarA_c.csv");
416  if(!inf){
417  if(world.rank() == 0) std::cerr << "Unable to open PbarA_c.csv" << std::endl;
418  exit(1);
419  }
420  getline(inf, strInput);
421  while(inf){
422  c.push_back(std::stod(strInput));
423  getline(inf, strInput);
424  }
425  inf.close();
426  //if(world.rank()==0) print("Done Reading PbarA_c");
427 
428  int n = c.size();
429  /*
430  for(int i=0; i < n; i++){
431  if(c[i]*exp(-t[i]*1e-30) < opthresh){
432  c.erase(c.begin()+i);
433  t.erase(t.begin()+i);
434  i--;
435  n--;
436  }
437  }
438  */
439  double fac = pow(2.0*pi,1.5);
440  for(int i = 0; i < n; i++){
441  c[i] = c[i]*fac;
442  }
443 
444  const Tensor<double> width = FunctionDefaults<3>::get_cell_width();
445  const int k = FunctionDefaults<3>::get_k();
446  double hi = width.normf(); // Diagonal width of cell
447 
448  std::vector<real_convolution_3d_ptr> gradG(3);
449 
450  for (int dir=0; dir<3; dir++) {
451  std::vector< ConvolutionND<double,3> > ops(n);
452  for (int mu=0; mu<n; mu++) {
453  // We cache the normalized operator so the factor is the value we must multiply
454  // by to recover the coeff we want.
455  double cc = std::pow(t[mu]/pi,1.5); // Normalization coeff
456  ops[mu].setfac(c[mu]/cc/width[dir]);
457 
458  for (int d=0; d<3; d++) {
459  if (d != dir) ops[mu].setop(d,GaussianConvolution1DCache<double>::get(k, t[mu]*width[d]*width[d], 0, false));
460  }
461  ops[mu].setop(dir,GaussianConvolution1DCache<double>::get(k, t[mu]*width[dir]*width[dir], 1, false));
462  }
463  gradG[dir] = real_convolution_3d_ptr(new SeparatedConvolution<double,3>(world, ops));
464  }
465 
466  return gradG;
467 
468 }
469 
470 #endif
Tbar
real_convolution_3d Tbar(World &world)
Definition: DKops.h:170
A2
real_convolution_3d A2(World &world)
Definition: DKops.h:286
Ebar
real_convolution_3d Ebar(World &world, double eps)
Definition: DKops.h:33
w
double w(double t, double eps)
Definition: DKops.h:22
Pbar
real_convolution_3d Pbar(World &world)
Definition: DKops.h:111
gradPbarA
std::vector< real_convolution_3d_ptr > gradPbarA(World &world)
Definition: DKops.h:397
PbarA
real_convolution_3d PbarA(World &world)
Definition: DKops.h:342
opthresh
static const double opthresh
Definition: DKops.h:16
q
double q(double t)
Definition: DKops.h:18
A
real_convolution_3d A(World &world)
Definition: DKops.h:230
Ebark
double Ebark(double t, double eps)
Definition: DKops.h:28
madness::BaseTensor::size
long size() const
Returns the number of elements in the tensor.
Definition: basetensor.h:138
madness::FunctionDefaults::get_k
static int get_k()
Returns the default wavelet order.
Definition: funcdefaults.h:266
madness::FunctionDefaults::get_cell_width
static const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition: funcdefaults.h:468
madness::Tensor::normf
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition: tensor.h:1726
madness::World
A parallel world class.
Definition: world.h:132
madness::World::rank
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition: world.h:318
constants.h
Defines common mathematical and physical constants.
R
static const double R
Definition: csqrt.cc:46
pow
static double pow(const double *a, const double *b)
Definition: lda.h:74
mra.h
Main include file for MADNESS and defines Function interface.
madness::constants::pi
const double pi
Mathematical constant .
Definition: constants.h:48
madness
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
madness::real_convolution_3d_ptr
std::shared_ptr< real_convolution_3d > real_convolution_3d_ptr
Definition: functypedefs.h:135
madness::real_convolution_3d
SeparatedConvolution< double, 3 > real_convolution_3d
Definition: functypedefs.h:121
madness::vec
Vector< T, sizeof...(Ts)+1 > vec(T t, Ts... ts)
Factory function for creating a madness::Vector.
Definition: vector.h:711
madness::dir
static const string dir
Definition: corepotential.cc:249
cc
const double cc
Definition: navstokes_cosines.cc:107
mu
const double mu
Definition: navstokes_cosines.cc:95
operator.h
Implements most functionality of separated operators.
c
static const double c
Definition: relops.cc:10
PI
static const double PI
Definition: relops.cc:12
k
static const long k
Definition: rk.cc:44
madness::GaussianConvolution1DCache
Definition: convolution1d.h:849
d
void d()
Definition: test_sig.cc:79
e
void e()
Definition: test_sig.cc:75
pi
static const double pi
Definition: testcosine.cc:6
Generated on Sat Oct 26 2024 16:26:20 for MADNESS by doxygen 1.9.1