1 /*
2 This file is part of MADNESS.
3
4 Copyright (C) 2007,2010 Oak Ridge National Laboratory
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20 For more information please contact:
21
22 Robert J. Harrison
23 Oak Ridge National Laboratory
24 One Bethel Valley Road
25 P.O. Box 2008, MS-6367
26
27 email: harrisonrj@ornl.gov
28 tel: 865-241-3937
29 fax: 865-572-0680
30
31
32 $Id$
33 */
34
35 /// \file testbsh.cc
36 /// \brief test the bsh operator
37
38 #include <madness/mra/mra.h>
39 #include <madness/mra/operator.h>
40 #include <madness/constants.h>
41
42 using namespace madness;
43
44 template <typename T, std::size_t NDIM>
45 class Gaussian : public FunctionFunctorInterface<T,NDIM> {
46 public:
47 typedef Vector<double,NDIM> coordT;
48 const coordT center;
49 const double exponent;
50 const T coefficient;
51
Gaussian(const coordT & center,double exponent,T coefficient)52 Gaussian(const coordT& center, double exponent, T coefficient)
53 : center(center), exponent(exponent), coefficient(coefficient) {};
54
operator ()(const coordT & x) const55 T operator()(const coordT& x) const {
56 double sum = 0.0;
57 for (std::size_t i=0; i<NDIM; ++i) {
58 double xx = center[i]-x[i];
59 sum += xx*xx;
60 };
61 return coefficient*exp(-exponent*sum);
62 };
63 };
64
65
66 double aa;
67
q(double r)68 double q(double r) {
69 double val;
70 if (r < 0.1e-4)
71 val = (0.2e1 * exp(0.1e1 / aa / 0.4e1) * exp(-0.1e1 / aa / 0.4e1) * sqrt(aa) / sqrt(constants::pi) + exp(0.1e1 / aa / 0.4e1) * erf(0.1e1 / sqrt(aa) / 0.2e1) - exp(0.1e1 / aa / 0.4e1) + (0.2e1 / 0.3e1 * exp(0.1e1 / aa / 0.4e1) * exp(-0.1e1 / aa / 0.4e1) * (0.1e1 / 0.2e1 - aa) * sqrt(aa) / sqrt(constants::pi) + exp(0.1e1 / aa / 0.4e1) * erf(0.1e1 / sqrt(aa) / 0.2e1) / 0.6e1 - exp(0.1e1 / aa / 0.4e1) / 0.6e1) * r * r);
72 else
73 val = ((-exp((0.1e1 + 0.4e1 * aa * r) / aa / 0.4e1) + exp(-(-0.1e1 + 0.4e1 * aa * r) / aa / 0.4e1) + exp((0.1e1 + 0.4e1 * aa * r) / aa / 0.4e1) * erf((0.2e1 * aa * r + 0.1e1) / sqrt(aa) / 0.2e1) + exp(-(-0.1e1 + 0.4e1 * aa * r) / aa / 0.4e1) * erf((-0.1e1 + 0.2e1 * aa * r) / sqrt(aa) / 0.2e1)) / r / 0.2e1);
74
75 return val / (4.0*constants::pi);
76 }
77
78
79 /// the result of the convolution
80 struct Qfunc : public FunctionFunctorInterface<double,3> {
operator ()Qfunc81 double operator()(const Vector<double,3>& x) const {
82 double r = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
83 return q(r);
84 }
85 };
86
87 template <typename T>
test_bsh(World & world)88 int test_bsh(World& world) {
89 double mu = 1.0;
90 std::vector<long> npt(3,201l);
91 typedef Vector<double,3> coordT;
92 typedef std::shared_ptr< FunctionFunctorInterface<T,3> > functorT;
93
94 int success=0;
95 if (world.rank() == 0)
96 print("Test BSH operation, type =",
97 archive::get_type_name<T>(),", ndim =",3);
98
99 FunctionDefaults<3>::set_cubic_cell(-100,100);
100 FunctionDefaults<3>::set_thresh(1e-6);
101 FunctionDefaults<3>::set_refine(true);
102 FunctionDefaults<3>::set_autorefine(true);
103 FunctionDefaults<3>::set_truncate_mode(1);
104 FunctionDefaults<3>::set_truncate_on_project(true);
105
106 for (int k=4; k<=12; k++) {
107 print("\n Testing with k", k);
108 FunctionDefaults<3>::set_k(k);
109
110 FunctionDefaults<3>::set_initial_level(6);
111 if (k >= 6) FunctionDefaults<3>::set_initial_level(5);
112
113 const coordT origin(0.0);
114 const double expnt = 100.0;
115 aa = expnt;
116 const double coeff = pow(expnt/constants::pi,1.5);
117
118 // the input function to be convolved
119 Function<T,3> f = FunctionFactory<T,3>(world).functor(functorT(new Gaussian<T,3>(origin, expnt, coeff)));
120 f.truncate();
121
122 double norm = f.trace();
123 double ferr = f.err(Gaussian<T,3>(origin, expnt, coeff));
124 if (world.rank() == 0) print("norm and error of the initial function", norm, ferr);
125
126 SeparatedConvolution<T,3> op = BSHOperator<3>(world, mu, 1e-4, 1e-8);
127 std::cout.precision(8);
128
129 // apply the convolution operator on the input function f
130 Function<T,3> ff = copy(f);
131 if (world.rank() == 0) print("applying - 1");
132 double start = cpu_time();
133 Function<T,3> opf = op(ff);
134 if (world.rank() == 0) print("done in time",cpu_time()-start);
135 ff.clear();
136 opf.verify_tree();
137 double opferr = opf.err(Qfunc());
138 if (world.rank() == 0) print("err in opf", opferr);
139 if (world.rank() == 0) print("err in f", ferr);
140
141 // here we are testing bsh, not the initial projection
142 //if ((opferr>ferr) and (opferr>FunctionDefaults<3>::get_thresh())) success++;
143 if (opferr>2*ferr) success++;
144
145 // //opf.truncate();
146 // Function<T,3> opinvopf = opf*(mu*mu);
147 // for (int axis=0; axis<3; ++axis) {
148 // opinvopf.gaxpy(1.0,diff(diff(opf,axis),axis).compress(),-1.0);
149 // }
150 // print("norm of (-del^2+mu^2)*G*f", opinvopf.norm2());
151 // Function<T,3> error = (f-opinvopf);
152 // print("error",error.norm2());
153
154 // opf.clear();
155 // opinvopf.clear();
156
157 // Function<T,3> g = (mu*mu)*f;
158 // for (int axis=0; axis<3; ++axis) {
159 // //g = g - diff(diff(f,axis),axis);
160 // }
161 // g = op(g);
162 // double derror=(g-f).norm2();
163 // print("norm of G*(-del^2+mu^2)*f",g.norm2());
164 // print("error",derror);
165 // if (derror> FunctionDefaults<3>::get_thresh()) success++;
166
167 }
168
169 world.gop.fence();
170 return success;
171 }
172
173
main(int argc,char ** argv)174 int main(int argc, char**argv) {
175 initialize(argc,argv);
176 World world(SafeMPI::COMM_WORLD);
177
178 int success=0;
179 try {
180 startup(world,argc,argv);
181
182 success=test_bsh<double>(world);
183
184 }
185 catch (const SafeMPI::Exception& e) {
186 print(e);
187 error("caught an MPI exception");
188 }
189 catch (const madness::MadnessException& e) {
190 print(e);
191 error("caught a MADNESS exception");
192 }
193 catch (const madness::TensorException& e) {
194 print(e);
195 error("caught a Tensor exception");
196 }
197 catch (char* s) {
198 print(s);
199 error("caught a c-string exception");
200 }
201 catch (const char* s) {
202 print(s);
203 error("caught a c-string exception");
204 }
205 catch (const std::string& s) {
206 print(s);
207 error("caught a string (class) exception");
208 }
209 catch (const std::exception& e) {
210 print(e.what());
211 error("caught an STL exception");
212 }
213 catch (...) {
214 error("caught unhandled exception");
215 }
216
217 world.gop.fence();
218 finalize();
219
220 return success;
221 }
222
223