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>
38 #include <madness/world/MADworld.h>
39 #include <madness/world/worlddc.h>
40 #include <madness/world/print.h>
41 #include <madness/misc/misc.h>
42 
43 #include <madness/tensor/tensor.h>
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 
54 namespace 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 
66 namespace 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:
271  typedef Tensor<T> tensorT ;
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 
284  bool is_second;
285  bool is_third;
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  }
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)
391  }
392  else {
393  if (bc_right == BC_DIRICHLET)
394  bdry_t.scale( pow(2.0,lev));
395  else if (bc_right ==BC_NEUMANN)
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(),
594  : DerivativeBase<T, NDIM>(world, axis, k, bc)
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 
612  void set_bspline1() {
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 
619  void set_bspline2() {
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 
626  void set_bspline3() {
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) {
692  set_is_second();
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>
706  }
707 
708 
709  /// Conveinence function returning derivative operator with periodic boundary conditions
710  template <typename T, std::size_t NDIM>
711  Derivative<T,NDIM>
714  }
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
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
const std::vector< long > vk
(k,...) used to initialize Tensors
Definition: derivative.h:80
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
Key< NDIM > neighbor(const keyT &key, int step) const
Definition: derivative.h:229
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
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
Future< argT > find_neighbor(const implT *f, const Key< NDIM > &key, int step) const
Definition: derivative.h:241
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:941
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
const std::shared_ptr< FunctionImpl< T, NDIM > > & get_impl() const
Returns a shared-pointer to the implementation.
Definition: mra.h:614
const Function< T, NDIM > & reconstruct(bool fence=true) const
Reconstructs the function, transforming into scaling function basis. Possible non-blocking comm.
Definition: mra.h:775
void set_impl(const std::shared_ptr< FunctionImpl< T, NDIM > > &impl)
Replace current FunctionImpl with provided new one.
Definition: mra.h:621
A future is a possibly yet unevaluated value.
Definition: future.h:373
remote_refT remote_ref(World &world) const
Returns a structure used to pass references to another process.
Definition: future.h:675
T & get(bool dowork=true) &
Gets the value, waiting if necessary.
Definition: future.h:574
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
static Key< NDIM > invalid()
Returns an invalid key.
Definition: key.h:104
Level level() const
Definition: key.h:159
const Vector< Translation, NDIM > & translation() const
Definition: key.h:164
bool is_invalid() const
Checks if a key is invalid.
Definition: key.h:109
static TaskAttributes hipri()
Definition: thread.h:450
A tensor is a multidimension array.
Definition: tensor.h:317
Tensor< T > cycledim(long nshift, long start, long end)
Returns new view/tensor cycling the sub-dimensions (start,...,end) with shift steps.
Definition: tensor.h:1641
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type scale(Q x)
Inplace multiplication by scalar of supported type (legacy name)
Definition: tensor.h:686
Makes a distributed container with specified attributes.
Definition: worlddc.h:866
ProcessID owner(const keyT &key) const
Returns processor that logically owns key (no communication)
Definition: worlddc.h:1034
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.
const double m
Definition: gfit.cc:199
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:190
#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:210
Header to declare stuff which has not yet found a home.
static const bool VERIFY_TREE
Definition: mra.h:57
File holds all helper structures necessary for the CC_Operator and CC2 class.
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
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
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
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d >>> &op, response_space &f)
Definition: basic_operators.cc:39
GenTensor< TENSOR_RESULT_TYPE(R, Q)> transform_dir(const GenTensor< R > &t, const Tensor< Q > &c, const int axis)
Definition: lowranktensor.h:1099
@ reconstructed
s coeffs at the leaves only
Definition: funcdefaults.h:59
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
response_space transpose(response_space &f)
Definition: basic_operators.cc:10
int64_t Translation
Definition: key.h:54
int Level
Definition: key.h:55
std::string get_mra_data_dir()
Definition: startup.cc:208
NDIM & f
Definition: mra.h:2416
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
double inner(response_space &a, response_space &b)
Definition: response_functions.h:442
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
Defines simple templates for printing to std::cout "a la Python".
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.
void d()
Definition: test_sig.cc:79
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