MADNESS 0.10.1
derivative.h
Go to the documentation of this file.
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#ifndef MADNESS_DERIVATIVE_H__INCLUDED
33#define MADNESS_DERIVATIVE_H__INCLUDED
34
35#include <iostream>
36#include <iomanip>
37#include <fstream>
40#include <madness/world/print.h>
41#include <madness/misc/misc.h>
42
45
46#include <madness/mra/key.h>
48
49
50/// \file mra/derivative.h
51/// \brief Declaration and initialization of tree traversal functions and generic derivative
52/// \ingroup mra
53
54namespace madness {
55
56 template<typename T, std::size_t NDIM>
57 class FunctionNode;
58
59 template<typename T, std::size_t NDIM>
60 class Function;
61
62}
63
64
65
66namespace madness {
67
68
69/// Tri-diagonal operator traversing tree primarily for derivative operator
70
71 /// \ingroup mra
72 template <typename T, std::size_t NDIM>
73 class DerivativeBase : public WorldObject< DerivativeBase<T, NDIM> > {
75 protected:
77 const std::size_t axis ; ///< Axis along which the operation is performed
78 const int k ; ///< Number of wavelets of the function
80 const std::vector<long> vk; ///< (k,...) used to initialize Tensors
81
82 public:
83 friend class FunctionImpl<T, NDIM>;
84
85 typedef Tensor<T> tensorT ; ///< regular tensors, like rm, etc
86 typedef GenTensor<T> coeffT ; ///< holding the node's coeffs (possibly low rank)
87 typedef Key<NDIM> keyT ;
88 typedef std::pair<keyT,coeffT> argT ;
93
94
97 , world(world)
98 , axis(axis)
99 , k(k)
100 , bc(bc)
101 , vk(NDIM,k)
102 {
103 // No! Cannot process incoming messages until the *derived* class is constructed.
104 // this->process_pending();
105 }
106
107 virtual ~DerivativeBase() { }
108
109 void forward_do_diff1(const implT* f, implT* df, const keyT& key,
110 const argT& left,
111 const argT& center,
112 const argT& right) const {
113
114 const dcT& coeffs = f->get_coeffs();
115 ProcessID owner = coeffs.owner(key);
116
117 if (owner == world.rank()) {
118 if (!left.second.has_data()) {
120 f, df, key, find_neighbor(f, key,-1), center, right,
122 }
123 else if (!right.second.has_data()) {
125 f, df, key, left, center, find_neighbor(f, key,1),
127 }
128 // Boundary node
129 else if (left.first.is_invalid() || right.first.is_invalid()) {
131 f, df, key, left, center, right);
132 }
133 // Interior node
134 else {
136 f, df, key, left, center, right);
137 }
138 }
139 else {
141 this, f, key, left, center, right, TaskAttributes::hipri());
142 }
143 }
144
145 void do_diff1(const implT* f, implT* df, const keyT& key,
146 const argT& left,
147 const argT& center,
148 const argT& right) const {
150
151// if (left.second.size()==0 || right.second.size()==0) {
152 if ((!left.second.has_data()) || (!right.second.has_data())) {
153 // One of the neighbors is below us in the tree ... recur down
154 df->get_coeffs().replace(key,nodeT(coeffT(),true));
155 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
156 const keyT& child = kit.key();
157 if ((child.translation()[axis]&1) == 0) {
158 // leftmost child automatically has right sibling
159 forward_do_diff1(f, df, child, left, center, center);
160 }
161 else {
162 // rightmost child automatically has left sibling
163 forward_do_diff1(f, df, child, center, center, right);
164 }
165 }
166 }
167 else {
168 forward_do_diff1(f, df, key, left, center, right);
169 }
170 }
171
172 virtual void do_diff2b(const implT* f, implT* df, const keyT& key,
173 const argT& left,
174 const argT& center,
175 const argT& right) const = 0;
176
177 virtual void do_diff2i(const implT* f, implT* df, const keyT& key,
178 const argT& left,
179 const argT& center,
180 const argT& right) const = 0;
181
182
183 /// Differentiate w.r.t. given coordinate (x=0, y=1, ...) with optional fence
184
185 /// Returns a new function with the same distribution
187 operator()(const functionT& f, bool fence=true) const {
188 if (VERIFY_TREE) f.verify_tree();
189 if (fence) f.change_tree_state(reconstructed);
190 MADNESS_CHECK_THROW(f.is_reconstructed(),"diff: trying to diff a compressed function without fencing");
191
192 functionT df;
193 df.set_impl(f,false);
194
195 df.get_impl()->diff(this, f.get_impl().get(), fence);
196 return df;
197 }
198
199
200 static bool enforce_bc(int bc_left, int bc_right, Level n, Translation& l) {
201 Translation two2n = 1ul << n;
202 if (l < 0) {
203 if (bc_left == BC_ZERO || bc_left == BC_FREE || bc_left == BC_DIRICHLET || bc_left == BC_ZERONEUMANN || bc_left == BC_NEUMANN) {
204 return false; // f=0 BC, or no BC, or nonzero f BC, or zero deriv BC, or nonzero deriv BC
205 }
206 else if (bc_left == BC_PERIODIC) {
207 l += two2n; // Periodic BC
208 MADNESS_ASSERT(bc_left == bc_right); //check that both BCs are periodic
209 }
210 else {
211 MADNESS_EXCEPTION("enforce_bc: confused left BC?",bc_left);
212 }
213 }
214 else if (l >= two2n) {
215 if (bc_right == BC_ZERO || bc_right == BC_FREE || bc_right == BC_DIRICHLET || bc_right == BC_ZERONEUMANN || bc_right == BC_NEUMANN) {
216 return false; // f=0 BC, or no BC, or nonzero f BC, or zero deriv BC, or nonzero deriv BC
217 }
218 else if (bc_right == BC_PERIODIC) {
219 l -= two2n; // Periodic BC
220 MADNESS_ASSERT(bc_left == bc_right); //check that both BCs are periodic
221 }
222 else {
223 MADNESS_EXCEPTION("enforce_bc: confused BC right?",bc_right);
224 }
225 }
226 return true;
227 }
228
229 Key<NDIM> neighbor(const keyT& key, int step) const {
231 l[axis] += step;
232 if (!enforce_bc(bc(axis,0), bc(axis,1), key.level(), l[axis])) {
233 return keyT::invalid();
234 }
235 else {
236 return keyT(key.level(),l);
237 }
238 }
239
241 find_neighbor(const implT* f, const Key<NDIM>& key, int step) const {
242 keyT neigh = neighbor(key, step);
243 if (neigh.is_invalid()) {
244 return Future<argT>(argT(neigh,coeffT(vk,f->get_tensor_args()))); // Zero bc
245 }
246 else {
247 Future<argT> result;
248 if (f->get_coeffs().is_local(neigh))
249 f->send(f->get_coeffs().owner(neigh), &implT::sock_it_to_me, neigh, result.remote_ref(world));
250 else
251 f->task(f->get_coeffs().owner(neigh), &implT::sock_it_to_me, neigh, result.remote_ref(world), TaskAttributes::hipri());
252 return result;
253 }
254 }
255
256
257 template <typename Archive> void serialize(const Archive& ar) const {
258 throw "NOT IMPLEMENTED";
259 }
260
261 }; // End of the DerivativeBase class
262
263
264 /// Implements derivatives operators with variety of boundary conditions on simulation domain
265 template <typename T, std::size_t NDIM>
266 class Derivative : public DerivativeBase<T, NDIM> {
267 private:
269
270 public:
272 typedef GenTensor<T> coeffT ; ///< holding the node's coeffs (possibly low rank)
273 typedef Key<NDIM> keyT ;
274 typedef std::pair<keyT,coeffT> argT ;
279
280 private:
281 const functionT g1; ///< Function describing the boundary condition on the right side
282 const functionT g2; ///< Function describing the boundary condition on the left side
283
286
287 // Tensors for holding the modified coefficients
288 Tensor<double> rm, r0, rp ; ///< Blocks of the derivative operator
289 Tensor<double> rmt, r0t, rpt ; ///< Blocks of the derivative operator, transposed
290 Tensor<double> left_rm, left_r0 ; ///< Blocks of the derivative for the left boundary
291 Tensor<double> left_rmt, left_r0t ; ///< Blocks of the derivative for the left boundary
292 Tensor<double> right_r0, right_rp; ///< Blocks of the derivative for the right boundary
293 Tensor<double> right_r0t, right_rpt; ///< Blocks of the derivative for the right boundary
294 Tensor<double> bv_left, bv_right ; ///< Blocks of the derivative operator for the boundary contribution
295
296
297 // Tensors for the bspline smoothed central difference operator
304
305 void do_diff2b(const implT* f, implT* df, const keyT& key,
306 const argT& left,
307 const argT& center,
308 const argT& right) const {
310 double lev = (double) key.level();
311
312 coeffT d;
313
314 //left boundary
315 if (l[this->axis] == 0) {
316
317 coeffT tensor_right=df->parent_to_child(right.second, right.first, this->neighbor(key,1));
318 coeffT tensor_center=df->parent_to_child(center.second, center.first, key);
319
320 d= transform_dir(tensor_right,left_rmt,this->axis);
321 d+=transform_dir(tensor_center,left_r0t,this->axis);
322 }
323 else {
324
325 coeffT tensor_left=df->parent_to_child(left.second, left.first, this->neighbor(key,-1));
326 coeffT tensor_center=df->parent_to_child(center.second, center.first, key);
327
328 d= transform_dir(tensor_left,right_rpt,this->axis);
329 d+=transform_dir(tensor_center,right_r0t,this->axis);
330 }
331
332 double fac = FunctionDefaults<NDIM>::get_rcell_width()[this->axis]*pow(2.0,lev);
333 if (is_second) fac *= fac;
334 else if (is_third) fac *= fac*fac;
335
336 d.scale(fac);
337 d.reduce_rank(df->get_thresh());
338 df->get_coeffs().replace(key,nodeT(d,false));
339
340
341 // This is the boundary contribution (formally in BoundaryDerivative)
342 int bc_left = this->bc(this->axis,0);
343 int bc_right = this->bc(this->axis,1);
344
345 Future<argT> found_argT;
346 tensorT bf, bdry_t;
347 //left boundary
348 if (l[this->axis] == 0) {
349 if (bc_left != BC_PERIODIC && bc_left != BC_FREE && bc_left != BC_ZERO && bc_left != BC_ZERONEUMANN) {
350 bf = copy(bv_left);
351 found_argT = g1.get_impl()->find_me(key);
352 }
353 else {
354 return;
355 }
356 }
357 else { //right boundary
358 if (bc_right != BC_PERIODIC && bc_right != BC_FREE && bc_right != BC_ZERO && bc_right != BC_ZERONEUMANN) {
359 bf = copy(bv_right);
360 found_argT = g2.get_impl()->find_me(key);
361 }
362 else {
363 return;
364 }
365 }
366#ifdef HAVE_PARSEC
367 std::cerr << "FATAL ERROR: PaRSEC does not support recursive task execution but Derivative::do_diff2b requires this. Use a different backend" << std::endl;
368 abort();
369#endif
370 const auto& found_argT_value = found_argT.get(); // do not recursively execute tasks to avoid making PaRSEC sad
371 tensorT gcoeffs = df->parent_to_child(found_argT_value.second, found_argT_value.first,key).full_tensor_copy();
372
373 //if (this->bc.get_bc().dim(0) == 1) {
374 if (NDIM == 1) {
375 bdry_t = gcoeffs[0]*bf;
376 }
377 else {
378 tensorT slice_aid(this->k); //vector of zeros
379 slice_aid[0] = 1;
380 tensorT tmp = inner(slice_aid, gcoeffs, 0, this->axis);
381 bdry_t = outer(bf,tmp);
382 if (this->axis) bdry_t = copy(bdry_t.cycledim(this->axis,0,this->axis)); // make it contiguous
383 }
384 bdry_t.scale(FunctionDefaults<NDIM>::get_rcell_width()[this->axis]);
385
386 if (l[this->axis]==0) {
387 if (bc_left == BC_DIRICHLET)
388 bdry_t.scale( pow(2.0,lev));
389 else if (bc_left ==BC_NEUMANN)
390 bdry_t.scale(FunctionDefaults<NDIM>::get_cell_width()[this->axis]);
391 }
392 else {
393 if (bc_right == BC_DIRICHLET)
394 bdry_t.scale( pow(2.0,lev));
395 else if (bc_right ==BC_NEUMANN)
396 bdry_t.scale(FunctionDefaults<NDIM>::get_cell_width()[this->axis]);
397 }
398
399 bdry_t += d.full_tensor_copy();;
400 df->get_coeffs().replace(key,nodeT(coeffT(bdry_t,df->get_thresh(),df->get_tensor_type()),false));
401 }
402
403 void do_diff2i(const implT* f, implT*df, const keyT& key,
404 const argT& left,
405 const argT& center,
406 const argT& right) const
407 {
408//#if !HAVE_GENTENSOR
409// coeffT d = madness::inner(rp,
410// df->parent_to_child(left.second, left.first, baseT::neighbor(key,-1)).swapdim(this->axis,0),
411// 1, 0);
412// inner_result(r0,
413// df->parent_to_child(center.second, center.first, key).swapdim(this->axis,0),
414// 1, 0, d);
415// inner_result(rm,
416// df->parent_to_child(right.second, right.first, baseT::neighbor(key,1)).swapdim(this->axis,0),
417// 1, 0, d);
418// // flo thinks this is wrong for higher dimensions -- need to cycledim
419// if (this->axis) d = copy(d.swapdim(this->axis,0)); // make it contiguous
420// d.scale(FunctionDefaults<NDIM>::get_rcell_width()[this->axis]*pow(2.0,(double) key.level()));
421// df->get_coeffs().replace(key,nodeT(d,false));
422//
423//#else
424 coeffT tensor_left=df->parent_to_child(left.second, left.first, this->neighbor(key,-1));
425 coeffT tensor_center=df->parent_to_child(center.second, center.first, key);
426 coeffT tensor_right=df->parent_to_child(right.second, right.first, this->neighbor(key,1));
427
428 coeffT d= transform_dir(tensor_left,rpt,this->axis);
429 d+=transform_dir(tensor_center,r0t,this->axis);
430 d+=transform_dir(tensor_right,rmt,this->axis);
431
432 double fac = FunctionDefaults<NDIM>::get_rcell_width()[this->axis]*pow(2.0,(double) key.level());
433 if (is_second) fac *= fac;
434 else if (is_third) fac *= fac*fac;
435
436 d.scale(fac);
437 d.reduce_rank(df->get_thresh());
438 df->get_coeffs().replace(key,nodeT(d,false));
439
440//#endif
441
442 }
443
445 is_second = false;
446 is_third = false;
447
448 r0 = Tensor<double>(this->k,this->k);
449 rp = Tensor<double>(this->k,this->k);
450 rm = Tensor<double>(this->k,this->k);
451
452 left_rm = Tensor<double>(this->k,this->k);
453 left_r0 = Tensor<double>(this->k,this->k);
454
455 right_r0 = Tensor<double>(this->k,this->k);
456 right_rp = Tensor<double>(this->k,this->k);
457
458 // These are the coefficients for the boundary contribution
459 bv_left = Tensor<double>(this->k);
460 bv_right = Tensor<double>(this->k);
461
462 int bc_left = this->bc(this->axis,0);
463 int bc_right = this->bc(this->axis,1);
464
465 double kphase = -1.0;
466 if (this->k%2 == 0) kphase = 1.0;
467 double iphase = 1.0;
468 for (int i=0; i<this->k; ++i) {
469 double jphase = 1.0;
470 for (int j=0; j<this->k; ++j) {
471 double gammaij = sqrt(double((2*i+1)*(2*j+1)));
472 double Kij;
473 if (((i-j)>0) && (((i-j)%2)==1))
474 Kij = 2.0;
475 else
476 Kij = 0.0;
477
478 r0(i,j) = 0.5*(1.0 - iphase*jphase - 2.0*Kij)*gammaij;
479 rm(i,j) = 0.5*jphase*gammaij;
480 rp(i,j) =-0.5*iphase*gammaij;
481
482 // Constraints on the derivative
483 if (bc_left == BC_ZERONEUMANN || bc_left == BC_NEUMANN) {
484 left_rm(i,j) = jphase*gammaij*0.5*(1.0 + iphase*kphase/this->k);
485
486 double phi_tmpj_left = 0;
487
488 for (int l=0; l<this->k; ++l) {
489 double gammalj = sqrt(double((2*l+1)*(2*j+1)));
490 double Klj;
491
492 if (((l-j)>0) && (((l-j)%2)==1)) Klj = 2.0;
493 else Klj = 0.0;
494
495 phi_tmpj_left += sqrt(double(2*l+1))*Klj*gammalj;
496 }
497 phi_tmpj_left = -jphase*phi_tmpj_left;
498 left_r0(i,j) = (0.5*(1.0 + iphase*kphase/this->k) - Kij)*gammaij + iphase*sqrt(double(2*i+1))*phi_tmpj_left/pow(this->k,2.);
499 }
500 else if (bc_left == BC_ZERO || bc_left == BC_DIRICHLET || bc_left == BC_FREE) {
501 left_rm(i,j) = rm(i,j);
502
503 // B.C. with a function
504 if (bc_left == BC_ZERO || bc_left == BC_DIRICHLET)
505 left_r0(i,j) = (0.5 - Kij)*gammaij;
506
507 // No B.C.
508 else if (bc_left == BC_FREE)
509 left_r0(i,j) = (0.5 - iphase*jphase - Kij)*gammaij;
510 }
511
512 // Constraints on the derivative
513 if (bc_right == BC_ZERONEUMANN || bc_right == BC_NEUMANN) {
514 right_rp(i,j) = -0.5*(iphase + kphase / this->k)*gammaij;
515
516 double phi_tmpj_right = 0;
517 for (int l=0; l<this->k; ++l) {
518 double gammalj = sqrt(double((2*l+1)*(2*j+1)));
519 double Klj;
520 if (((l-j)>0) && (((l-j)%2)==1)) Klj = 2.0;
521 else Klj = 0.0;
522 phi_tmpj_right += sqrt(double(2*l+1))*Klj*gammalj;
523 }
524 right_r0(i,j) = -(0.5*jphase*(iphase+ kphase/this->k) + Kij)*gammaij + sqrt(double(2*i+1))*phi_tmpj_right/pow(this->k,2.);
525 }
526 else if (bc_right == BC_ZERO || bc_right == BC_FREE || bc_right == BC_DIRICHLET) {
527 right_rp(i,j) = rp(i,j);
528
529 // Zero BC
530 if (bc_right == BC_ZERO || bc_right == BC_DIRICHLET)
531 right_r0(i,j) = -(0.5*iphase*jphase + Kij)*gammaij;
532
533 // No BC
534 else if (bc_right == BC_FREE)
535 right_r0(i,j) = (1.0 - 0.5*iphase*jphase - Kij)*gammaij;
536
537 }
538
539 jphase = -jphase;
540 }
541 iphase = -iphase;
542 }
543
544 // Coefficients for the boundary contributions
545 iphase = 1.0;
546 for (int i=0; i<this->k; ++i) {
547 iphase = -iphase;
548
549 if (bc_left == BC_DIRICHLET)
550 bv_left(i) = iphase*sqrt(double(2*i+1)); // vector for left dirichlet BC
551 else if(bc_left == BC_NEUMANN)
552 bv_left(i) = -iphase*sqrt(double(2*i+1))/pow(this->k,2.); // vector for left deriv BC
553 else
554 bv_left(i) = 0.0;
555
556 if (bc_right == BC_DIRICHLET)
557 bv_right(i) = sqrt(double(2*i+1)); // vector for right dirichlet BC
558 else if (bc_right == BC_NEUMANN)
559 bv_right(i) = sqrt(double(2*i+1))/pow(this->k,2.); // vector for right deriv BC
560 else
561 bv_right(i) = 0.0;
562 }
563
564 r0t = transpose(r0);
565 rpt = transpose(rp);
566 rmt = transpose(rm);
567
570
573
574 //print(rm.normf(),r0.normf(),rp.normf(),left_rm.normf(),left_r0.normf(),right_r0.normf(),right_rp.normf(),bv_left.normf(),bv_right.normf());
575 }
576
577 public:
578 typedef T opT;
579
580 /// Constructs a derivative operator
581
582 /// @param world The world
583 /// @param axis The direction to differentiate
584 /// @param bc Boundary conditions (default from FunctionDefaults)
585 /// @param g1 Function providing left boundary value (default empty)
586 /// @param g2 Function providing right boundary value (default empty)
587 /// @param k Wavelet order (default from FunctionDefaults)
589 std::size_t axis,
591 const functionT g1=functionT(),
592 const functionT g2=functionT(),
595 , g1(g1)
596 , g2(g2)
597 {
600 g1.reconstruct();
601 g2.reconstruct();
602
603 this->process_pending();
604 }
605
606 virtual ~Derivative() { }
607
608 void set_is_first() {is_second = false; is_third = false;}
609 void set_is_second() {is_second = true; is_third=false;}
610 void set_is_third() {is_second = false; is_third = true;}
611
614 if(k > 18) throw "Bspline derivatives are only available up to k=18";
615 std::string filename = get_mra_data_dir() + "/b-spline-deriv1.txt";
617 }
618
621 if(k > 18) throw "Bspline derivatives are only available up to k=18";
622 std::string filename = get_mra_data_dir() + "/b-spline-deriv2.txt";
624 }
625
628 if(k > 18) throw "Bspline derivatives are only available up to k=18";
629 std::string filename = get_mra_data_dir() + "/b-spline-deriv3.txt";
631 }
632
633 void set_ble1() {
635 if(k > 15) throw "BLE derivatives are only available up to k=15";
636 std::string filename = get_mra_data_dir() + "/ble-first.txt";
638 }
639
640 void set_ble2() {
642 if(k > 15) throw "BLE derivatives are only available up to k=15";
643 std::string filename = get_mra_data_dir() + "/ble-second.txt";
645 }
646
647 void read_from_file(const std::string& filename, unsigned int order = 1) {
648
649 Tensor<double> r0_bsp(this->k,this->k);
650 Tensor<double> rp_bsp(this->k,this->k);
651 Tensor<double> rm_bsp(this->k,this->k);
652
653 std::ifstream f(filename);
654 bool found=false;
655
656 for (int m; f >> m; ) {
657 if (m == this->k) {
658 for (int i=0; i<m; i++)
659 for (int j=0; j<m; j++)
660 MADNESS_CHECK(f >> rp_bsp(i,j));
661 for (int i=0; i<m; i++)
662 for (int j=0; j<m; j++)
663 MADNESS_CHECK(f >> r0_bsp(i,j));
664 for (int i=0; i<m; i++)
665 for (int j=0; j<m; j++)
666 MADNESS_CHECK(f >> rm_bsp(i,j));
667 found = true;
668 break;
669 }
670 else {
671 double junk;
672 for (int i=0; i<3*m*m; i++)
673 MADNESS_CHECK(f >> junk);
674 }
675 }
676 MADNESS_CHECK(found);
680
682
684
686
687 // Get scaling factor right for higher order derivatives
688 if (order == 1) {
689 set_is_first();
690 }
691 else if(order == 2) {
693 }
694 else if(order == 3) {
695 set_is_third();
696 }
697 }
698 };
699
700
701 /// Convenience function returning derivative operator with free-space boundary conditions
702 template <typename T, std::size_t NDIM>
703 Derivative<T,NDIM>
707
708
709 /// Conveinence function returning derivative operator with periodic boundary conditions
710 template <typename T, std::size_t NDIM>
711 Derivative<T,NDIM>
715
716 /// Applies derivative operator to function (for syntactic equivalence to integral operator apply)
717 template <typename T, std::size_t NDIM>
718 Function<T,NDIM>
719 apply(const Derivative<T,NDIM>& D, const Function<T,NDIM>& f, bool fence=true) {
720 return D(f,fence);
721 }
722
723 /// Convenience function returning vector of derivative operators implementing grad (\f$ \nabla \f$)
724
725 /// This will only work for BC_ZERO, BC_PERIODIC, BC_FREE and
726 /// BC_ZERONEUMANN since we are not passing in any boundary
727 /// functions.
728 template <typename T, std::size_t NDIM>
729 std::vector< std::shared_ptr< Derivative<T,NDIM> > >
733 std::vector< std::shared_ptr< Derivative<T,NDIM> > > r(NDIM);
734 for (std::size_t d=0; d<NDIM; ++d) {
735 MADNESS_CHECK(bc(d,0)!=BC_DIRICHLET && bc(d,1)!=BC_DIRICHLET);
736 MADNESS_CHECK(bc(d,0)!=BC_NEUMANN && bc(d,1)!=BC_NEUMANN);
737 r[d].reset(new Derivative<T,NDIM>(world,d,bc,Function<T,NDIM>(),Function<T,NDIM>(),k));
738 }
739 return r;
740 }
741
742
743 namespace archive {
744 template <class Archive, class T, std::size_t NDIM>
745 struct ArchiveLoadImpl<Archive,const DerivativeBase<T,NDIM>*> {
746 static void load(const Archive& ar, const DerivativeBase<T,NDIM>*& ptr) {
748 ar & p;
749 ptr = static_cast< const DerivativeBase<T,NDIM>* >(p);
750 }
751 };
752
753 template <class Archive, class T, std::size_t NDIM>
754 struct ArchiveStoreImpl<Archive,const DerivativeBase<T,NDIM>*> {
755 static void store(const Archive& ar, const DerivativeBase<T,NDIM>* const & ptr) {
756 ar & ptr->id();
757 }
758 };
759 }
760
761} // End of the madness namespace
762
763#endif // MADNESS_MRA_DERIVATIVE_H_INCLUDED
This header should include pretty much everything needed for the parallel runtime.
This class is used to specify boundary conditions for all operators.
Definition funcdefaults.h:101
Tri-diagonal operator traversing tree primarily for derivative operator.
Definition derivative.h:73
void do_diff1(const implT *f, implT *df, const keyT &key, const argT &left, const argT &center, const argT &right) const
Definition derivative.h:145
GenTensor< T > coeffT
holding the node's coeffs (possibly low rank)
Definition derivative.h:86
static bool enforce_bc(int bc_left, int bc_right, Level n, Translation &l)
Definition derivative.h:200
DerivativeBase(World &world, std::size_t axis, int k, BoundaryConditions< NDIM > bc)
Definition derivative.h:95
Key< NDIM > keyT
Definition derivative.h:87
const BoundaryConditions< NDIM > bc
Definition derivative.h:79
Tensor< T > tensorT
regular tensors, like rm, etc
Definition derivative.h:85
const std::vector< long > vk
(k,...) used to initialize Tensors
Definition derivative.h:80
Key< NDIM > neighbor(const keyT &key, int step) const
Definition derivative.h:229
WorldContainer< Key< NDIM >, FunctionNode< T, NDIM > > dcT
Definition derivative.h:91
virtual ~DerivativeBase()
Definition derivative.h:107
FunctionImpl< T, NDIM > implT
Definition derivative.h:89
FunctionNode< T, NDIM > nodeT
Definition derivative.h:92
Function< T, NDIM > functionT
Definition derivative.h:90
void forward_do_diff1(const implT *f, implT *df, const keyT &key, const argT &left, const argT &center, const argT &right) const
Definition derivative.h:109
const int k
Number of wavelets of the function.
Definition derivative.h:78
WorldObject< DerivativeBase< T, NDIM > > woT
Definition derivative.h:74
void serialize(const Archive &ar) const
Definition derivative.h:257
Future< argT > find_neighbor(const implT *f, const Key< NDIM > &key, int step) const
Definition derivative.h:241
Function< T, NDIM > operator()(const functionT &f, bool fence=true) const
Differentiate w.r.t. given coordinate (x=0, y=1, ...) with optional fence.
Definition derivative.h:187
virtual void do_diff2i(const implT *f, implT *df, const keyT &key, const argT &left, const argT &center, const argT &right) const =0
const std::size_t axis
Axis along which the operation is performed.
Definition derivative.h:77
World & world
Definition derivative.h:76
virtual void do_diff2b(const implT *f, implT *df, const keyT &key, const argT &left, const argT &center, const argT &right) const =0
std::pair< keyT, coeffT > argT
Definition derivative.h:88
Implements derivatives operators with variety of boundary conditions on simulation domain.
Definition derivative.h:266
Tensor< double > right_r0t
Definition derivative.h:293
void set_ble2()
Definition derivative.h:640
Tensor< double > rmt
Definition derivative.h:289
Tensor< double > bv_left
Definition derivative.h:294
void set_bspline1()
Definition derivative.h:612
Tensor< double > r0
Definition derivative.h:288
Tensor< double > rp_bsp
Definition derivative.h:300
bool is_second
Definition derivative.h:284
Tensor< double > right_rp
Blocks of the derivative for the right boundary.
Definition derivative.h:292
void set_is_second()
Definition derivative.h:609
Derivative(World &world, std::size_t axis, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), const functionT g1=functionT(), const functionT g2=functionT(), int k=FunctionDefaults< NDIM >::get_k())
Constructs a derivative operator.
Definition derivative.h:588
Function< T, NDIM > functionT
Definition derivative.h:276
Tensor< double > r0t
Definition derivative.h:289
std::pair< keyT, coeffT > argT
Definition derivative.h:274
FunctionImpl< T, NDIM > implT
Definition derivative.h:275
Tensor< double > right_rpt
Blocks of the derivative for the right boundary.
Definition derivative.h:293
Tensor< double > left_rmt
Definition derivative.h:291
Tensor< double > rp_bsp_t
Definition derivative.h:303
virtual ~Derivative()
Definition derivative.h:606
GenTensor< T > coeffT
holding the node's coeffs (possibly low rank)
Definition derivative.h:272
const functionT g2
Function describing the boundary condition on the left side.
Definition derivative.h:282
void read_from_file(const std::string &filename, unsigned int order=1)
Definition derivative.h:647
Tensor< double > rp
Blocks of the derivative operator.
Definition derivative.h:288
void do_diff2i(const implT *f, implT *df, const keyT &key, const argT &left, const argT &center, const argT &right) const
Definition derivative.h:403
bool is_third
Definition derivative.h:285
void set_bspline3()
Definition derivative.h:626
Tensor< double > rm_bsp
Definition derivative.h:299
void set_bspline2()
Definition derivative.h:619
Tensor< double > rm
Definition derivative.h:288
void initCoefficients()
Definition derivative.h:444
Tensor< double > left_r0
Blocks of the derivative for the left boundary.
Definition derivative.h:290
Tensor< double > rpt
Blocks of the derivative operator, transposed.
Definition derivative.h:289
void do_diff2b(const implT *f, implT *df, const keyT &key, const argT &left, const argT &center, const argT &right) const
Definition derivative.h:305
void set_is_third()
Definition derivative.h:610
T opT
Definition derivative.h:578
Tensor< double > rm_bsp_t
Definition derivative.h:302
Tensor< double > left_r0t
Blocks of the derivative for the left boundary.
Definition derivative.h:291
Tensor< T > tensorT
Definition derivative.h:271
void set_is_first()
Definition derivative.h:608
FunctionNode< T, NDIM > nodeT
Definition derivative.h:278
Tensor< double > bv_right
Blocks of the derivative operator for the boundary contribution.
Definition derivative.h:294
Key< NDIM > keyT
Definition derivative.h:273
const functionT g1
Function describing the boundary condition on the right side.
Definition derivative.h:281
void set_ble1()
Definition derivative.h:633
Tensor< double > left_rm
Definition derivative.h:290
WorldContainer< Key< NDIM >, FunctionNode< T, NDIM > > dcT
Definition derivative.h:277
Tensor< double > r0_bsp
Definition derivative.h:298
Tensor< double > right_r0
Definition derivative.h:292
DerivativeBase< T, NDIM > baseT
Definition derivative.h:268
Tensor< double > r0_bsp_t
Definition derivative.h:301
FunctionDefaults holds default paramaters as static class members.
Definition funcdefaults.h:204
static int get_k()
Returns the default wavelet order.
Definition funcdefaults.h:266
static const Tensor< double > & get_rcell_width()
Returns the reciprocal of the width of each user cell dimension.
Definition funcdefaults.h:473
FunctionImpl holds all Function state to facilitate shallow copy semantics.
Definition funcimpl.h:942
void sock_it_to_me(const keyT &key, const RemoteReference< FutureImpl< std::pair< keyT, coeffT > > > &ref) const
Walk up the tree returning pair(key,node) for first node with coefficients.
Definition mraimpl.h:2847
double get_thresh() const
Definition mraimpl.h:307
TensorType get_tensor_type() const
Definition mraimpl.h:298
const coeffT parent_to_child(const coeffT &s, const keyT &parent, const keyT &child) const
Directly project parent coeffs to child coeffs.
Definition mraimpl.h:3178
const dcT & get_coeffs() const
Definition mraimpl.h:322
FunctionNode holds the coefficients, etc., at each node of the 2^NDIM-tree.
Definition funcimpl.h:124
A multiresolution adaptive numerical function.
Definition mra.h:122
A future is a possibly yet unevaluated value.
Definition future.h:373
T & get(bool dowork=true) &
Gets the value, waiting if necessary.
Definition future.h:574
remote_refT remote_ref(World &world) const
Returns a structure used to pass references to another process.
Definition future.h:675
Definition lowranktensor.h:59
Tensor< T > full_tensor_copy() const
Definition gentensor.h:206
Iterates in lexical order thru all children of a key.
Definition key.h:374
Key is the index for a node of the 2^NDIM-tree.
Definition key.h:66
Level level() const
Definition key.h:159
bool is_invalid() const
Checks if a key is invalid.
Definition key.h:109
const Vector< Translation, NDIM > & translation() const
Definition key.h:164
static Key< NDIM > invalid()
Returns an invalid key.
Definition key.h:104
static TaskAttributes hipri()
Definition thread.h:450
A tensor is a multidimension array.
Definition tensor.h:317
A simple, fixed dimension vector.
Definition vector.h:64
Makes a distributed container with specified attributes.
Definition worlddc.h:866
void replace(const pairT &datum)
Inserts/replaces key+value pair (non-blocking communication if key not local)
Definition worlddc.h:974
Implements most parts of a globally addressable object (via unique ID).
Definition world_object.h:364
const uniqueidT & id() const
Returns the globally unique object ID.
Definition world_object.h:711
void process_pending()
To be called from derived constructor to process pending messages.
Definition world_object.h:656
detail::task_result_type< memfnT >::futureT task(ProcessID dest, memfnT memfn, const TaskAttributes &attr=TaskAttributes()) const
Sends task to derived class method returnT (this->*memfn)().
Definition world_object.h:1005
A parallel world class.
Definition world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition world.h:318
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
Provides FunctionDefaults and utilities for coordinate transformation.
auto T(World &world, response_space &f) -> response_space
Definition global_functions.cc:34
Multidimension Key for MRA tree and associated iterators.
static double pow(const double *a, const double *b)
Definition lda.h:74
#define MADNESS_CHECK(condition)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:182
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition madness_exception.h:119
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition madness_exception.h:134
#define MADNESS_CHECK_THROW(condition, msg)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:207
Header to declare stuff which has not yet found a home.
static const bool VERIFY_TREE
Definition mra.h:57
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
@ BC_DIRICHLET
Definition funcdefaults.h:56
@ BC_NEUMANN
Definition funcdefaults.h:56
@ BC_ZERO
Definition funcdefaults.h:56
@ BC_PERIODIC
Definition funcdefaults.h:56
@ BC_ZERONEUMANN
Definition funcdefaults.h:56
@ BC_FREE
Definition funcdefaults.h:56
static const char * filename
Definition legendre.cc:96
@ reconstructed
s coeffs at the leaves only
Definition funcdefaults.h:59
response_space transpose(response_space &f)
Definition basic_operators.cc:10
int64_t Translation
Definition key.h:54
std::vector< std::shared_ptr< Derivative< T, NDIM > > > gradient_operator(World &world, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int k=FunctionDefaults< NDIM >::get_k())
Convenience function returning vector of derivative operators implementing grad ( )
Definition derivative.h:730
Derivative< T, NDIM > periodic_derivative(World &world, int axis, int k=FunctionDefaults< NDIM >::get_k())
Conveinence function returning derivative operator with periodic boundary conditions.
Definition derivative.h:712
int Level
Definition key.h:55
std::enable_if< std::is_base_of< ProjectorBase, projT >::value, OuterProjector< projT, projQ > >::type outer(const projT &p0, const projQ &p1)
Definition projector.h:457
std::string get_mra_data_dir()
Definition startup.cc:208
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d > > > &op, response_space &f)
Definition basic_operators.cc:39
Derivative< T, NDIM > free_space_derivative(World &world, int axis, int k=FunctionDefaults< NDIM >::get_k())
Convenience function returning derivative operator with free-space boundary conditions.
Definition derivative.h:704
NDIM & f
Definition mra.h:2416
double inner(response_space &a, response_space &b)
Definition response_functions.h:442
GenTensor< TENSOR_RESULT_TYPE(R, Q)> transform_dir(const GenTensor< R > &t, const Tensor< Q > &c, const int axis)
Definition lowranktensor.h:1099
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition mra.h:2002
static const double d
Definition nonlinschro.cc:121
Defines simple templates for printing to std::cout "a la Python".
static const double m
Definition relops.cc:9
static const long k
Definition rk.cc:44
Definition test_ar.cc:204
static void load(const Archive &ar, const DerivativeBase< T, NDIM > *&ptr)
Definition derivative.h:746
Default load of an object via serialize(ar, t).
Definition archive.h:666
static void store(const Archive &ar, const DerivativeBase< T, NDIM > *const &ptr)
Definition derivative.h:755
Default store of an object via serialize(ar, t).
Definition archive.h:611
Defines and implements most of Tensor.
static const std::size_t NDIM
Definition testpdiff.cc:42
std::size_t axis
Definition testpdiff.cc:59
Implements WorldContainer.
int ProcessID
Used to clearly identify process number/rank.
Definition worldtypes.h:43