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