32 #ifndef MADNESS_DERIVATIVE_H__INCLUDED
33 #define MADNESS_DERIVATIVE_H__INCLUDED
56 template<
typename T, std::
size_t NDIM>
59 template<
typename T, std::
size_t NDIM>
72 template <
typename T, std::
size_t NDIM>
80 const std::vector<long>
vk;
88 typedef std::pair<keyT,coeffT>
argT ;
112 const argT& right)
const {
114 const dcT& coeffs =
f->get_coeffs();
118 if (!left.second.has_data()) {
123 else if (!right.second.has_data()) {
129 else if (left.first.is_invalid() || right.first.is_invalid()) {
131 f, df, key, left, center, right);
136 f, df, key, left, center, right);
148 const argT& right)
const {
152 if ((!left.second.has_data()) || (!right.second.has_data())) {
156 const keyT& child = kit.key();
175 const argT& right)
const = 0;
180 const argT& right)
const = 0;
190 MADNESS_CHECK_THROW(
f.is_reconstructed(),
"diff: trying to diff a compressed function without fencing");
195 df.
get_impl()->diff(
this,
f.get_impl().get(), fence);
214 else if (l >= two2n) {
248 if (
f->get_coeffs().is_local(neigh))
257 template <
typename Archive>
void serialize(
const Archive& ar)
const {
258 throw "NOT IMPLEMENTED";
265 template <
typename T, std::
size_t NDIM>
274 typedef std::pair<keyT,coeffT>
argT ;
308 const argT& right)
const {
310 double lev = (double) key.
level();
315 if (l[this->
axis] == 0) {
317 coeffT tensor_right=df->
parent_to_child(right.second, right.first, this->neighbor(key,1));
325 coeffT tensor_left=df->
parent_to_child(left.second, left.first, this->neighbor(key,-1));
342 int bc_left = this->
bc(this->
axis,0);
343 int bc_right = this->
bc(this->
axis,1);
348 if (l[this->
axis] == 0) {
367 std::cerr <<
"FATAL ERROR: PaRSEC does not support recursive task execution but Derivative::do_diff2b requires this. Use a different backend" << std::endl;
370 const auto& found_argT_value = found_argT.
get();
375 bdry_t = gcoeffs[0]*bf;
381 bdry_t =
outer(bf,tmp);
386 if (l[this->
axis]==0) {
399 bdry_t +=
d.full_tensor_copy();;
406 const argT& right)
const
462 int bc_left = this->
bc(this->
axis,0);
463 int bc_right = this->
bc(this->
axis,1);
465 double kphase = -1.0;
466 if (this->
k%2 == 0) kphase = 1.0;
468 for (
int i=0; i<this->
k; ++i) {
470 for (
int j=0; j<this->
k; ++j) {
471 double gammaij = sqrt(
double((2*i+1)*(2*j+1)));
473 if (((i-j)>0) && (((i-j)%2)==1))
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;
484 left_rm(i,j) = jphase*gammaij*0.5*(1.0 + iphase*kphase/this->
k);
486 double phi_tmpj_left = 0;
488 for (
int l=0; l<this->
k; ++l) {
489 double gammalj = sqrt(
double((2*l+1)*(2*j+1)));
492 if (((l-j)>0) && (((l-j)%2)==1)) Klj = 2.0;
495 phi_tmpj_left += sqrt(
double(2*l+1))*Klj*gammalj;
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.);
505 left_r0(i,j) = (0.5 - Kij)*gammaij;
509 left_r0(i,j) = (0.5 - iphase*jphase - Kij)*gammaij;
514 right_rp(i,j) = -0.5*(iphase + kphase / this->
k)*gammaij;
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)));
520 if (((l-j)>0) && (((l-j)%2)==1)) Klj = 2.0;
522 phi_tmpj_right += sqrt(
double(2*l+1))*Klj*gammalj;
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.);
531 right_r0(i,j) = -(0.5*iphase*jphase + Kij)*gammaij;
535 right_r0(i,j) = (1.0 - 0.5*iphase*jphase - Kij)*gammaij;
546 for (
int i=0; i<this->
k; ++i) {
550 bv_left(i) = iphase*sqrt(
double(2*i+1));
552 bv_left(i) = -iphase*sqrt(
double(2*i+1))/
pow(this->k,2.);
614 if(
k > 18)
throw "Bspline derivatives are only available up to k=18";
621 if(
k > 18)
throw "Bspline derivatives are only available up to k=18";
628 if(
k > 18)
throw "Bspline derivatives are only available up to k=18";
635 if(
k > 15)
throw "BLE derivatives are only available up to k=15";
642 if(
k > 15)
throw "BLE derivatives are only available up to k=15";
656 for (
int m;
f >>
m; ) {
658 for (
int i=0; i<
m; i++)
659 for (
int j=0; j<
m; j++)
661 for (
int i=0; i<
m; i++)
662 for (
int j=0; j<
m; j++)
664 for (
int i=0; i<
m; i++)
665 for (
int j=0; j<
m; j++)
672 for (
int i=0; i<3*
m*
m; i++)
691 else if(order == 2) {
694 else if(order == 3) {
702 template <
typename T, std::
size_t NDIM>
710 template <
typename T, std::
size_t NDIM>
717 template <
typename T, std::
size_t NDIM>
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) {
744 template <
class Archive,
class T, std::
size_t NDIM>
753 template <
class Archive,
class T, std::
size_t NDIM>
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 ¢er, 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 ¢er, 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 ¢er, 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 ¢er, 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 ¢er, 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 ¢er, 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