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/tensor/solvers.h>
41 #include <madness/constants.h>
42
43 using namespace madness;
44
45 typedef std::shared_ptr< WorldDCPmapInterface< Key<3> > > pmapT;
46 typedef Vector<double,3> coordT;
47 typedef std::shared_ptr< FunctionFunctorInterface<double,3> > functorT;
48 typedef Function<double,3> functionT;
49 typedef std::vector<functionT> vecfuncT;
50 typedef std::pair<vecfuncT,vecfuncT> pairvecfuncT;
51 typedef std::vector<pairvecfuncT> subspaceT;
52 typedef Tensor<double> tensorT;
53 typedef FunctionFactory<double,3> factoryT;
54 typedef SeparatedConvolution<double,3> operatorT;
55 typedef std::shared_ptr<operatorT> poperatorT;
56
57 // The delta fn approx is (1/sqrt(2*pi*s^2))*exp(-x^2 / (2*s^2))
58
59 class Sphere : public FunctionFunctorInterface<double,3> {
60 public:
61 const coordT center;
62 const double radius;
63 const double sigma;
64 const double expnt;
65 const double norm;
66
Sphere(const coordT & center,double radius,double sigma)67 Sphere(const coordT& center, double radius, double sigma)
68 : center(center)
69 , radius(radius)
70 , sigma(sigma)
71 , expnt(0.5/(sigma*sigma))
72 , norm(1.0/(sigma*sqrt(2.0*constants::pi)))
73 {}
74
operator ()(const coordT & x) const75 double operator()(const coordT& x) const {
76 double xx = (x[0]-center[0]);
77 double yy = (x[1]-center[1]);
78 double zz = (x[2]-center[2]);
79 double r = sqrt(xx*xx + yy*yy + zz*zz);
80 double dist = r - radius; // Normal distance from boundary
81 double arg = dist*dist*expnt; // Argument to exponential
82 if (arg > 35.0)
83 return 0.0;
84 else
85 return norm*exp(-arg);
86 }
87 };
88
89
90 class DSphere : public FunctionFunctorInterface<double,3> {
91 public:
92 const coordT center;
93 const double radius;
94 const double sigma;
95 const double expnt;
96 const double norm;
97
DSphere(const coordT & center,double radius,double sigma)98 DSphere(const coordT& center, double radius, double sigma)
99 : center(center)
100 , radius(radius)
101 , sigma(sigma)
102 , expnt(0.5/(sigma*sigma))
103 , norm(1.0/(sigma*sqrt(2.0*constants::pi)))
104 {}
105
operator ()(const coordT & x) const106 double operator()(const coordT& x) const {
107 double xx = (x[0]-center[0]);
108 double yy = (x[1]-center[1]);
109 double zz = (x[2]-center[2]);
110 double r = sqrt(xx*xx + yy*yy + zz*zz);
111 double dist = r - radius; // Normal distance from boundary
112 double arg = dist*dist*expnt; // Argument to exponential
113 if (arg > 35.0)
114 return 0.0;
115 else
116 return -2.0*dist*expnt*norm*exp(-arg);
117 }
118 };
119
120
main(int argc,char ** argv)121 int main(int argc, char**argv) {
122 initialize(argc,argv);
123 World world(SafeMPI::COMM_WORLD);
124 try {
125 startup(world,argc,argv);
126 std::cout.precision(6);
127
128 // MADNESS simulation parameters
129 const int k = 6;
130 const double thresh = 1e-5;
131 const double L = 4.0;
132
133 // Boundary width
134 double sigma = 0.1;
135
136 // Plotting info
137 int npt = 3001;
138 coordT lo(0.0), hi(0.0);
139 lo[0] = -L;
140 hi[0] = L;
141
142 FunctionDefaults<3>::set_cubic_cell(-L,L);
143 FunctionDefaults<3>::set_k(k);
144 FunctionDefaults<3>::set_thresh(thresh);
145 FunctionDefaults<3>::set_initial_level(3);
146 FunctionDefaults<3>::set_refine(true);
147 FunctionDefaults<3>::set_autorefine(false);
148 FunctionDefaults<3>::set_truncate_mode(1);
149 FunctionDefaults<3>::set_truncate_on_project(true);
150
151 // Inner and outer shells
152 functionT spha = factoryT(world).functor(functorT(new Sphere(coordT(0.0), 1.0, sigma)));
153 functionT sphb = factoryT(world).functor(functorT(new Sphere(coordT(0.0), 3.0, sigma)));
154
155 //functionT dspha = factoryT(world).functor(functorT(new DSphere(coordT(0.0), 1.0, sigma)));
156 //functionT dsphb = factoryT(world).functor(functorT(new DSphere(coordT(0.0), 3.0, sigma)));
157
158 print("inner shell norm", spha.norm2());
159 print("outer shell norm", sphb.norm2());
160
161 functionT phi = 3.0*spha + sphb; // Delta function at the boundary
162 functionT u0 = -3.0*spha + sphb; // The desired boundary conditions * phi
163 //functionT du0= -1.0*dspha + dsphb; // The desired boundary conditions * phi
164 //dspha.clear(); dsphb.clear();
165 spha.clear(); sphb.clear();
166
167 phi.truncate(thresh*0.1);
168 u0.truncate(thresh*0.1);
169 //du0.truncate(thresh*0.1);
170
171 plot_line("phi.dat", npt, lo, hi, phi);
172 plot_line("u0.dat", npt, lo, hi, u0);
173
174 operatorT op = CoulombOperator(world, sigma*0.1, thresh);
175 // Starting values
176 double mu = 0.25;
177 functionT f = u0; //
178 //du0.clear();
179 // functionT lambda = f; // lambda is actually lambdaphi
180
181 functionT u; // This will be current solution
182 functionT c; // This will be the value of the constraint
183
184 for (int muiter=0; muiter<5; ++muiter) {
185 // for (int lamiter=0; lamiter<20; ++lamiter) {
186 vecfuncT rvec;
187 vecfuncT fvec;
188 std::vector<double> rnorms;
189 tensorT Q;
190 for (int m=0; m<20; ++m) {
191 f.truncate();
192 fvec.push_back(f);
193 u = op(fvec[m]);
194 u.scale(0.25/constants::pi);
195 u.truncate(); // Current soln
196
197 c = u*phi - u0; // Constraint violation
198 double cnorm = c.norm2();
199
200 // rvec.push_back(fvec[m] - lambda + (phi*c)*(1.0/mu)); // Residual AugLag
201 rvec.push_back(fvec[m] + phi*c*(1.0/mu)); // Residual QuadPen
202 rvec.back().truncate();
203 double rnorm = rvec[m].norm2();
204 rnorms.push_back(rnorm);
205
206 // Update subspace matrix
207 tensorT Qnew(m+1,m+1);
208 if (m>0) {
209 Qnew(Slice(0,m-1),Slice(0,m-1)) = Q;
210 }
211 Q = Qnew;
212 for (int i=0; i<=m; ++i) {
213 Q(i,m) = fvec[i].inner(rvec[m]);
214 if (m != i) Q(m,i) = fvec[m].inner(rvec[i]);
215 }
216 // print("Subspace matrix");
217 // print(Q);
218
219 tensorT coeff = KAIN(Q);
220
221 // print("Subspace solution");
222 // print(coeff);
223
224 f = coeff[0]*(fvec[0]-rvec[0]);
225 for (int i=1; i<=m; ++i) {
226 f.gaxpy(1.0,fvec[i]-rvec[i],coeff[i]);
227 }
228
229 // Step restriction
230 double snorm = (f - fvec[m]).norm2();
231 double fnorm = fvec[m].norm2();
232
233 if (snorm > fnorm*0.3) {
234 double damp = fnorm*0.3/snorm;
235 print(" damping", damp);
236 f.gaxpy(damp, fvec[m], 1.0-damp);
237 }
238 f.truncate();
239
240 print(" iteration", m, "mu", mu, "fnorm", fnorm, "snorm", snorm, "rnorm", rnorm, "cnorm", cnorm);
241
242 if (snorm < 0.1*cnorm) {
243 // Current solution will be in f
244 break;
245 }
246 }
247
248 char fname[256];
249 int lamiter = 0;
250 sprintf(fname, "c-%2.2d-%2.2d.dat", lamiter, muiter);
251 plot_line(fname, npt, lo, hi, c);
252 sprintf(fname, "u-%2.2d-%2.2d.dat", lamiter, muiter);
253 plot_line(fname, npt, lo, hi, u);
254 sprintf(fname, "f-%2.2d-%2.2d.dat", lamiter, muiter);
255 plot_line(fname, npt, lo, hi, f);
256 // sprintf(fname, "l-%2.2d-%2.2d.dat", lamiter, muiter);
257 // plot_line(fname, npt, lo, hi, lambda);
258
259 // double damp = 1.0;
260 // functionT lambdanew = lambda - phi* c * (damp/mu);
261
262 // print("\nupdating lambda", lamiter, (lambdanew-lambda).norm2());
263 // lambda = lambdanew;
264 // lambda.truncate();
265 // }
266 mu *= 0.5;
267 // lambda = f + phi * c * (1.0/mu);
268 print("\nupdating mu", muiter, mu);
269 }
270
271 }
272 catch (const SafeMPI::Exception& e) {
273 print(e);
274 error("caught an MPI exception");
275 }
276 catch (const madness::MadnessException& e) {
277 print(e);
278 error("caught a MADNESS exception");
279 }
280 catch (const madness::TensorException& e) {
281 print(e);
282 error("caught a Tensor exception");
283 }
284 catch (char* s) {
285 print(s);
286 error("caught a c-string exception");
287 }
288 catch (const char* s) {
289 print(s);
290 error("caught a c-string exception");
291 }
292 catch (const std::string& s) {
293 print(s);
294 error("caught a string (class) exception");
295 }
296 catch (const std::exception& e) {
297 print(e.what());
298 error("caught an STL exception");
299 }
300 catch (...) {
301 error("caught unhandled exception");
302 }
303
304 world.gop.fence();
305 finalize();
306
307 return 0;
308 }
309
310