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_MRA_FUNCDEFAULTS_H__INCLUDED 33 #define MADNESS_MRA_FUNCDEFAULTS_H__INCLUDED 34 35 36 /// \file funcdefaults.h 37 /// \brief Provides FunctionDefaults and utilities for coordinate transformation 38 /// \ingroup mrabcext 39 40 #include <madness/world/MADworld.h> 41 #include <madness/world/vector.h> 42 #include <madness/world/worlddc.h> 43 #include <madness/tensor/tensor.h> 44 #include <madness/mra/key.h> 45 46 namespace madness { 47 template <typename T, std::size_t NDIM> class FunctionImpl; 48 49 /// The maximum wavelet order presently supported 50 static const int MAXK = 30; 51 52 /// The maximum depth of refinement possible 53 static const int MAXLEVEL = 8*sizeof(Translation)-2; 54 55 enum BCType {BC_ZERO, BC_PERIODIC, BC_FREE, BC_DIRICHLET, BC_ZERONEUMANN, BC_NEUMANN}; 56 57 /*! 58 \brief This class is used to specify boundary conditions for all operators 59 \ingroup mrabcext 60 61 Exterior boundary conditions (i.e., on the simulation domain) 62 are associated with operators (not functions). The types of 63 boundary conditions available are in the enum BCType. 64 65 The default boundary conditions are obtained from the FunctionDefaults. 66 For non-zero Dirichlet and Neumann conditions additional information 67 must be provided when derivative operators are constructed. For integral 68 operators, only periodic and free space are supported. 69 */ 70 template<std::size_t NDIM> 71 class BoundaryConditions { 72 private: 73 // Used to use STL vector but static data on a MAC was 74 // causing problems. 75 BCType bc[NDIM*2]; 76 77 public: 78 /// Constructor. Default boundary condition set to free space 79 BoundaryConditions(BCType code=BC_FREE) 80 { 81 for (std::size_t i=0; i<NDIM*2; ++i) bc[i] = code; 82 } 83 84 /// Copy constructor is deep BoundaryConditions(const BoundaryConditions<NDIM> & other)85 BoundaryConditions(const BoundaryConditions<NDIM>& other) 86 { 87 *this = other; 88 } 89 90 /// Assignment makes deep copy 91 BoundaryConditions<NDIM>& 92 operator=(const BoundaryConditions<NDIM>& other) { 93 if (&other != this) { 94 for (std::size_t i=0; i<NDIM*2; ++i) bc[i] = other.bc[i]; 95 } 96 return *this; 97 } 98 99 /// Returns value of boundary condition 100 101 /// @param d Dimension (0,...,NDIM-1) for boundary condition 102 /// @param i Side (0=left, 1=right) for boundary condition 103 /// @return Value of boundary condition operator()104 BCType operator()(std::size_t d, int i) const { 105 MADNESS_ASSERT(d<NDIM && i>=0 && i<2); 106 return bc[2*d+i]; 107 } 108 109 /// Returns reference to boundary condition 110 111 /// @param d Dimension (0,...,NDIM-1) for boundary condition 112 /// @param i Side (0=left, 1=right) for boundary condition 113 /// @return Value of boundary condition operator()114 BCType& operator()(std::size_t d, int i) { 115 MADNESS_ASSERT(d<NDIM && i>=0 && i<2); 116 return bc[2*d+i]; 117 } 118 119 template <typename Archive> serialize(const Archive & ar)120 void serialize(const Archive& ar) { 121 ar & bc; 122 } 123 124 /// Translates code into human readable string 125 126 /// @param code Code for boundary condition 127 /// @return String describing boundary condition code code_as_string(BCType code)128 static const char* code_as_string(BCType code) { 129 static const char* codes[] = {"zero","periodic","free","Dirichlet","zero Neumann","Neumann"}; 130 return codes[code]; 131 } 132 133 /// Convenience for application of integral operators 134 135 /// @return Returns a vector indicating if each dimension is periodic is_periodic()136 std::vector<bool> is_periodic() const { 137 std::vector<bool> v(NDIM); 138 for (std::size_t d=0; d<NDIM; ++d) 139 v[d] = (bc[2*d]==BC_PERIODIC); 140 return v; 141 } 142 }; 143 144 145 146 template <std::size_t NDIM> 147 static 148 inline 149 std::ostream& operator << (std::ostream& s, const BoundaryConditions<NDIM>& bc) { 150 s << "BoundaryConditions("; 151 for (unsigned int d=0; d<NDIM; ++d) { 152 s << bc.code_as_string(bc(d,0)) << ":" << bc.code_as_string(bc(d,1)); 153 if (d == NDIM-1) 154 s << ")"; 155 else 156 s << ", "; 157 } 158 return s; 159 } 160 161 162 /// FunctionDefaults holds default paramaters as static class members 163 164 /// Declared and initialized in mra.cc and/or funcimpl::initialize. 165 /// 166 /// Currently all functions of the same dimension share the same cell dimensions 167 /// since they are stored inside FunctionDefaults ... if you change the 168 /// cell dimensions *all* functions of that dimension are affected. 169 /// 170 /// N.B. Ultimately, we may need to make these defaults specific to each 171 /// world, as should be all global state. 172 /// \ingroup mra 173 template <std::size_t NDIM> 174 class FunctionDefaults { 175 private: 176 static int k; ///< Wavelet order 177 static double thresh; ///< Truncation threshold 178 static int initial_level; ///< Initial level for fine scale projection 179 static int special_level; ///< Minimum level for fine scale projection of special boxes 180 static int max_refine_level; ///< Level at which to stop refinement 181 static int truncate_mode; ///< Truncation method 182 static bool refine; ///< Whether to refine new functions 183 static bool autorefine; ///< Whether to autorefine in multiplication, etc. 184 static bool debug; ///< Controls output of debug info 185 static bool truncate_on_project; ///< If true initial projection inserts at n-1 not n 186 static bool apply_randomize; ///< If true use randomization for load balancing in apply integral operator 187 static bool project_randomize; ///< If true use randomization for load balancing in project/refine 188 static BoundaryConditions<NDIM> bc; ///< Default boundary conditions 189 static Tensor<double> cell ; ///< cell[NDIM][2] Simulation cell, cell(0,0)=xlo, cell(0,1)=xhi, ... 190 static Tensor<double> cell_width;///< Width of simulation cell in each dimension 191 static Tensor<double> rcell_width; ///< Reciprocal of width 192 static double cell_volume; ///< Volume of simulation cell 193 static double cell_min_width; ///< Size of smallest dimension 194 static TensorType tt; ///< structure of the tensor in FunctionNode 195 static std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > > pmap; ///< Default mapping of keys to processes 196 recompute_cell_info()197 static void recompute_cell_info() { 198 MADNESS_ASSERT(cell.dim(0)==NDIM && cell.dim(1)==2 && cell.ndim()==2); 199 cell_width = cell(_,1)-cell(_,0); 200 cell_volume = cell_width.product(); 201 cell_min_width = cell_width.min(); 202 rcell_width = copy(cell_width); 203 for (std::size_t i=0; i<NDIM; ++i) rcell_width(i) = 1.0/rcell_width(i); 204 } 205 206 public: 207 208 MADNESS_PRAGMA_CLANG(diagnostic push) 209 MADNESS_PRAGMA_CLANG(diagnostic ignored "-Wundefined-var-template") 210 211 /// Used to set defaults to k=7, thresh=1-5, for a unit cube [0,1]. 212 static void set_defaults(World& world); 213 214 static void print(); 215 216 /// Returns the default wavelet order get_k()217 static int get_k() { 218 return k; 219 } 220 221 /// Sets the default wavelet order 222 223 /// Existing functions are unaffacted. set_k(int value)224 static void set_k(int value) { 225 k=value; 226 MADNESS_ASSERT(k>0 && k<=MAXK); 227 } 228 229 /// Returns the default threshold get_thresh()230 static const double& get_thresh() { 231 return thresh; 232 } 233 234 /// Sets the default threshold 235 236 /// Existing functions are unaffected set_thresh(double value)237 static void set_thresh(double value) { 238 thresh=value; 239 } 240 241 /// Returns the default initial projection level get_initial_level()242 static int get_initial_level() { 243 return initial_level; 244 } 245 246 /// Returns the default projection level for special boxes get_special_level()247 static int get_special_level() { 248 return special_level; 249 } 250 251 /// Sets the default initial projection level 252 253 /// Existing functions are unaffected set_initial_level(int value)254 static void set_initial_level(int value) { 255 initial_level=value; 256 MADNESS_ASSERT(value>0 && value<MAXLEVEL); 257 } 258 259 /// Existing functions are unaffected set_special_level(int value)260 static void set_special_level(int value) { 261 special_level=value; 262 MADNESS_ASSERT(value>=0 && value<MAXLEVEL); 263 MADNESS_ASSERT(max_refine_level>=special_level); 264 } 265 266 /// Gets the default maximum adaptive refinement level get_max_refine_level()267 static int get_max_refine_level() { 268 return max_refine_level; 269 } 270 271 /// Sets the default maximum adaptive refinement level 272 273 /// Existing functions are unaffected set_max_refine_level(int value)274 static void set_max_refine_level(int value) { 275 max_refine_level=value; 276 MADNESS_ASSERT(value>0 && value<MAXLEVEL); 277 MADNESS_ASSERT(max_refine_level>=initial_level); 278 MADNESS_ASSERT(max_refine_level>=special_level); 279 } 280 281 /// Gets the default truncation mode get_truncate_mode()282 static int get_truncate_mode() { 283 return truncate_mode; 284 } 285 286 /// Sets the default truncation mode 287 288 /// Existing functions are unaffected set_truncate_mode(int value)289 static void set_truncate_mode(int value) { 290 truncate_mode=value; 291 MADNESS_ASSERT(value>=0 && value<4); 292 } 293 294 /// Gets the default adaptive refinement flag get_refine()295 static bool get_refine() { 296 return refine; 297 } 298 299 /// Sets the default adaptive refinement flag 300 301 /// Existing functions are unaffected set_refine(bool value)302 static void set_refine(bool value) { 303 refine=value; 304 } 305 306 /// Gets the default adaptive autorefinement flag get_autorefine()307 static bool get_autorefine() { 308 return autorefine; 309 } 310 311 /// Sets the default adaptive autorefinement flag 312 313 /// Existing functions are unaffected set_autorefine(bool value)314 static void set_autorefine(bool value) { 315 autorefine=value; 316 } 317 318 /// Gets the default debug flag (is this used anymore?) get_debug()319 static bool get_debug() { 320 return debug; 321 } 322 323 /// Sets the default debug flag (is this used anymore?) 324 325 /// Not sure if this does anything useful set_debug(bool value)326 static void set_debug(bool value) { 327 debug=value; 328 } 329 330 /// Gets the default truncate on project flag get_truncate_on_project()331 static bool get_truncate_on_project() { 332 return truncate_on_project; 333 } 334 335 /// Sets the default truncate on project flag 336 337 /// Existing functions are unaffected set_truncate_on_project(bool value)338 static void set_truncate_on_project(bool value) { 339 truncate_on_project=value; 340 } 341 342 /// Gets the random load balancing for integral operators flag get_apply_randomize()343 static bool get_apply_randomize() { 344 return apply_randomize; 345 } 346 347 /// Sets the random load balancing for integral operators flag set_apply_randomize(bool value)348 static void set_apply_randomize(bool value) { 349 apply_randomize=value; 350 } 351 352 353 /// Gets the random load balancing for projection flag get_project_randomize()354 static bool get_project_randomize() { 355 return project_randomize; 356 } 357 358 /// Sets the random load balancing for projection flag set_project_randomize(bool value)359 static void set_project_randomize(bool value) { 360 project_randomize=value; 361 } 362 363 /// Returns the default boundary conditions get_bc()364 static const BoundaryConditions<NDIM>& get_bc() { 365 return bc; 366 } 367 368 /// Sets the default boundary conditions set_bc(const BoundaryConditions<NDIM> & value)369 static void set_bc(const BoundaryConditions<NDIM>& value) { 370 bc=value; 371 } 372 373 /// Returns the default tensor type get_tensor_type()374 static TensorType get_tensor_type() { 375 return tt; 376 } 377 378 /// Sets the default tensor type set_tensor_type(const TensorType & t)379 static void set_tensor_type(const TensorType& t) { 380 #if HAVE_GENTENSOR 381 tt=t; 382 #else 383 tt=TT_FULL; 384 #endif 385 } 386 387 /// adapt the special level to resolve the smallest length scale 388 static int set_length_scale(const double lo,const size_t k=get_k()) { 389 const double dk = (double) k; 390 double Lmax=FunctionDefaults<NDIM>::get_cell_width().max(); 391 double lo_sim=lo/Lmax; // lo in simulation coordinates; 392 const int special_level=Level(-log2(lo_sim*dk)); 393 return special_level; 394 } 395 396 /// Gets the user cell for the simulation get_cell()397 static const Tensor<double>& get_cell() { 398 return cell; 399 } 400 401 /// Gets the user cell for the simulation 402 403 /// Existing functions are probably rendered useless set_cell(const Tensor<double> & value)404 static void set_cell(const Tensor<double>& value) { 405 cell=copy(value); 406 recompute_cell_info(); 407 } 408 409 /// Sets the user cell to be cubic with each dimension having range \c [lo,hi] 410 411 /// Existing functions are probably rendered useless set_cubic_cell(double lo,double hi)412 static void set_cubic_cell(double lo, double hi) { 413 cell(_,0)=lo; 414 cell(_,1)=hi; 415 recompute_cell_info(); 416 } 417 418 /// Returns the width of each user cell dimension get_cell_width()419 static const Tensor<double>& get_cell_width() { 420 return cell_width; 421 } 422 423 /// Returns the reciprocal of the width of each user cell dimension get_rcell_width()424 static const Tensor<double>& get_rcell_width() { 425 return rcell_width; 426 } 427 428 /// Returns the minimum width of any user cell dimension get_cell_min_width()429 static double get_cell_min_width() { 430 return cell_min_width; 431 } 432 433 /// Returns the volume of the user cell get_cell_volume()434 static double get_cell_volume() { 435 return cell_volume; 436 } 437 438 /// Returns the default process map get_pmap()439 static std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > >& get_pmap() { 440 return pmap; 441 } 442 443 /// Sets the default process map (does \em not redistribute existing functions) 444 445 /// Existing functions are probably rendered useless set_pmap(const std::shared_ptr<WorldDCPmapInterface<Key<NDIM>>> & value)446 static void set_pmap(const std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > >& value) { 447 pmap = value; 448 } 449 450 /// Sets the default process map and redistributes all functions using the old map redistribute(World & world,const std::shared_ptr<WorldDCPmapInterface<Key<NDIM>>> & newpmap)451 static void redistribute(World& world, const std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > >& newpmap) { 452 pmap->redistribute(world,newpmap); 453 pmap = newpmap; 454 } 455 456 MADNESS_PRAGMA_CLANG(diagnostic pop) 457 458 }; 459 460 461 /// Convert user coords (cell[][]) to simulation coords ([0,1]^ndim) 462 template <std::size_t NDIM> user_to_sim(const Vector<double,NDIM> & xuser,Vector<double,NDIM> & xsim)463 static inline void user_to_sim(const Vector<double,NDIM>& xuser, Vector<double,NDIM>& xsim) { 464 for (std::size_t i=0; i<NDIM; ++i) 465 xsim[i] = (xuser[i] - FunctionDefaults<NDIM>::get_cell()(i,0)) * FunctionDefaults<NDIM>::get_rcell_width()[i]; 466 } 467 468 /// Returns the box at level n that contains the given point in simulation coordinates 469 /// @param[in] pt point in simulation coordinates 470 /// @param[in] n the level of the box 471 template <typename T, std::size_t NDIM> simpt2key(const Vector<T,NDIM> & pt,Level n)472 static inline Key<NDIM> simpt2key(const Vector<T,NDIM>& pt, Level n){ 473 Vector<Translation,NDIM> l; 474 double twon = std::pow(2.0, double(n)); 475 for (std::size_t i=0; i<NDIM; ++i) { 476 l[i] = Translation(twon*pt[i]); 477 } 478 return Key<NDIM>(n,l); 479 } 480 481 /// Convert simulation coords ([0,1]^ndim) to user coords (FunctionDefaults<NDIM>::get_cell()) 482 template <std::size_t NDIM> sim_to_user(const Vector<double,NDIM> & xsim,Vector<double,NDIM> & xuser)483 static void sim_to_user(const Vector<double,NDIM>& xsim, Vector<double,NDIM>& xuser) { 484 for (std::size_t i=0; i<NDIM; ++i) 485 xuser[i] = xsim[i]*FunctionDefaults<NDIM>::get_cell_width()[i] + FunctionDefaults<NDIM>::get_cell()(i,0); 486 } 487 488 489 } 490 #endif // MADNESS_MRA_FUNCDEFAULTS_H__INCLUDED 491