diff options
| -rw-r--r-- | Carpet/Carpet/README | 5 | ||||
| -rw-r--r-- | Carpet/Carpet/src/Shutdown.cc | 47 | ||||
| -rw-r--r-- | Carpet/Carpet/src/carpet_public.hh | 299 | ||||
| -rw-r--r-- | Carpet/Carpet/src/make.code.defn | 3 | ||||
| -rw-r--r-- | Carpet/CarpetIOASCII/README | 3 | ||||
| -rw-r--r-- | Carpet/CarpetInterp/README | 3 | ||||
| -rw-r--r-- | Carpet/CarpetLib/README | 3 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/schedule.ccl | 9 | ||||
| -rw-r--r-- | Carpet/CarpetSlab/README | 3 | ||||
| -rw-r--r-- | Carpet/CarpetSlab/src/slab.h | 22 | ||||
| -rw-r--r-- | Carpet/CarpetTest/src/carpettest_check_arguments.F77 | 16 | ||||
| -rw-r--r-- | CarpetAttic/CarpetIOFlexIO/README | 4 | ||||
| -rw-r--r-- | CarpetDev/CarpetCG/src/CG.cc | 872 | ||||
| -rw-r--r-- | CarpetDev/CarpetJacobi/param.ccl | 14 | ||||
| -rw-r--r-- | CarpetDev/CarpetJacobi/src/Jacobi.cc | 698 |
15 files changed, 867 insertions, 1134 deletions
diff --git a/Carpet/Carpet/README b/Carpet/Carpet/README index 9b38c07cd..81065fb42 100644 --- a/Carpet/Carpet/README +++ b/Carpet/Carpet/README @@ -1,8 +1,9 @@ Cactus Code Thorn Carpet Authors : Erik Schnetter <schnetter@uni-tuebingen.de> -CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/README,v 1.1 2001/03/01 13:40:10 eschnett Exp $ +CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/README,v 1.2 2004/01/25 14:57:27 schnetter Exp $ -------------------------------------------------------------------------- Purpose of the thorn: -This thorn provides a parallel AMR (adaptive mesh refinement) driver with MPI. +This thorn provides a parallel AMR (adaptive mesh refinement) driver +with MPI. diff --git a/Carpet/Carpet/src/Shutdown.cc b/Carpet/Carpet/src/Shutdown.cc index cc01fbd06..b4d79ac2d 100644 --- a/Carpet/Carpet/src/Shutdown.cc +++ b/Carpet/Carpet/src/Shutdown.cc @@ -10,7 +10,7 @@ #include "carpet.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Shutdown.cc,v 1.12 2003/08/03 17:09:02 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Shutdown.cc,v 1.13 2004/01/25 14:57:27 schnetter Exp $"; CCTK_FILEVERSION(Carpet_Carpet_Shutdown_cc); } @@ -26,40 +26,45 @@ namespace Carpet { { DECLARE_CCTK_PARAMETERS; - Waypoint ("starting Shutdown..."); + Waypoint ("Starting shutdown"); const int convlev = 0; cGH* cgh = fc->GH[convlev]; - Waypoint ("Current time is %g", cgh->cctk_time); - - BEGIN_REFLEVEL_LOOP(cgh) { - BEGIN_MGLEVEL_LOOP(cgh) { - - do_global_mode = reflevel == 0; + for (int rl=reflevels-1; rl>=0; --rl) { + BEGIN_REVERSE_MGLEVEL_LOOP(cgh) { + enter_level_mode (cgh, rl); + do_global_mode = reflevel==0; + do_meta_mode = do_global_mode && mglevel==mglevels-1; - Waypoint ("%*sCurrent time is %g%s", 2*reflevel, "", - cgh->cctk_time, - do_global_mode ? " (global time)" : ""); + Checkpoint ("Shutdown at iteration %d time %g%s%s", + cgh->cctk_iteration, (double)cgh->cctk_time, + (do_global_mode ? " (global)" : ""), + (do_meta_mode ? " (meta)" : "")); // Terminate - Waypoint ("%*sScheduling TERMINATE", 2*reflevel, ""); + Checkpoint ("Scheduling TERMINATE"); CCTK_ScheduleTraverse ("CCTK_TERMINATE", cgh, CallFunction); - - } END_MGLEVEL_LOOP; - } END_REFLEVEL_LOOP; - - do_global_mode = true; + + leave_level_mode (cgh); + } END_REVERSE_MGLEVEL_LOOP; + } // for rl - // Shutdown - Waypoint ("Scheduling SHUTDOWN"); - CCTK_ScheduleTraverse ("CCTK_SHUTDOWN", cgh, CallFunction); + BEGIN_REVERSE_MGLEVEL_LOOP(cgh) { + do_global_mode = true; + do_meta_mode = mglevel==mglevels-1; + + // Shutdown + Checkpoint ("Scheduling SHUTDOWN"); + CCTK_ScheduleTraverse ("CCTK_SHUTDOWN", cgh, CallFunction); + + } END_REVERSE_MGLEVEL_LOOP; CCTK_PRINTSEPARATOR; printf ("Done.\n"); // earlier checkpoint before finalising MPI - Waypoint ("done with Shutdown."); + Waypoint ("Done with shutdown"); dist::finalize(); diff --git a/Carpet/Carpet/src/carpet_public.hh b/Carpet/Carpet/src/carpet_public.hh index 64ce2703b..fbe9d8c0f 100644 --- a/Carpet/Carpet/src/carpet_public.hh +++ b/Carpet/Carpet/src/carpet_public.hh @@ -1,302 +1,15 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.39 2003/11/05 16:18:37 schnetter Exp $ - -// It is assumed that the number of components of all arrays is equal -// to the number of components of the grid functions, and that their -// distribution onto the processors is the same, and that all -// processors own the same number of components. +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.40 2004/01/25 14:57:27 schnetter Exp $ #ifndef CARPET_PUBLIC_HH #define CARPET_PUBLIC_HH -#include <vector> - -#include "cctk.h" -#include "cctk_Schedule.h" - -#include "dh.hh" -#include "gf.hh" -#include "ggf.hh" -#include "gh.hh" -#include "th.hh" - -#include "carpet_public.h" - - - // Stuff with C linkage #include "carpet_public.h" - - -namespace Carpet { - - - - const int dim = 3; - - - - // Handle from CCTK_RegisterGHExtension - extern int GHExtension; - - // Refinement factor - extern int reffact; - - // Refinement factor on finest grid - extern int maxreflevelfact; - - // Multigrid levels - extern int mglevels; - - // Multigrid factor - extern int mgfact; - - // Maximum number of refinement levels - extern int maxreflevels; - - // Multigrid factor on coarsest grid - extern int maxmglevelfact; - - - - // Current position on the grid hierarchy - extern int reflevel; - extern int mglevel; - extern int component; - - // Current refinement factor - extern int reflevelfact; - - // Current multigrid factor - extern int mglevelfact; - - // Is this the time for a global mode call? - extern bool do_global_mode; - - // Is prolongation enabled? - extern bool do_prolongate; - - // Current times on the refinement levels - extern vector<CCTK_REAL> refleveltimes; - extern CCTK_REAL delta_time; - - - - // Data for grid functions - - // The grid hierarchy - extern gh<dim>* hh; - extern th<dim>* tt; - extern dh<dim>* dd; - - // Data for everything - struct arrdesc { - // points to hh etc. for GF, and is unique for SCALAR and ARRAY - cGroupDynamicData info; - gh<dim>* hh; - th<dim>* tt; - dh<dim>* dd; - vector<ggf<dim>*> data; // [var] - bool do_transfer; // prolongate and restrict - // VGF - }; - extern vector<arrdesc> arrdata; // [group] - - - - // Checksums - struct ckdesc { - bool valid; - unsigned int sum; - }; - // [rl][c][group][var][tl] - extern vector<vector<vector<vector<vector<ckdesc> > > > > checksums; - - - - // Registered functions - void* SetupGH (tFleshConfig* fc, int convLevel, cGH* cgh); - - int Initialise (tFleshConfig* config); - int Evolve (tFleshConfig* config); - int Shutdown (tFleshConfig* config); - int CallFunction (void* function, cFunctionData* attribute, void* data); - - int SyncGroup (const cGH* cgh, const char* groupname); - int EnableGroupStorage (const cGH* cgh, const char* groupname); - int DisableGroupStorage (const cGH* cgh, const char* groupname); - int EnableGroupComm (const cGH* cgh, const char* groupname); - int DisableGroupComm (const cGH* cgh, const char* groupname); - int Barrier (const cGH* cgh); - int Exit (cGH* cgh, int retval); - int Abort (cGH* cgh, int retval); - int MyProc (const cGH* cgh); - int nProcs (const cGH* cgh); - const int* ArrayGroupSizeB (const cGH* cgh, int dir, int group, - const char* groupname); - int QueryGroupStorageB (const cGH* cgh, int group, const char* groupname); - int GroupDynamicData (const cGH* cgh, int group, cGroupDynamicData* data); - - - - // Functions for recomposing the grid hierarchy - void RegisterRegridRoutine (int (*routine)(const cGH * cckgGH, - gh<dim>::rexts& bbsss, - gh<dim>::rbnds& obss, - gh<dim>::rprocs& pss)); - - void SplitRegions (const cGH* cgh, vector<bbox<int,dim> >& bbs, - vector<vect<vect<bool,2>,dim> >& obs); - void SplitRegions_AlongZ (const cGH* cgh, vector<bbox<int,dim> >& bbs, - vector<vect<vect<bool,2>,dim> >& obs); - void SplitRegions_AlongDir (const cGH* cgh, vector<bbox<int,dim> >& bbs, - vector<vect<vect<bool,2>,dim> >& obs, - const int dir); - void SplitRegions_Automatic (const cGH* cgh, vector<bbox<int,dim> >& bbs, - vector<vect<vect<bool,2>,dim> >& obs); - - void MakeProcessors (const cGH* cgh, const gh<dim>::rexts& bbsss, - gh<dim>::rprocs& pss); - - - - // Helper functions - void set_reflevel (cGH* cgh, int rl); - void set_mglevel (cGH* cgh, int ml); - void set_component (cGH* cgh, int c); - - - - // Refinement level iterator - -#define BEGIN_REFLEVEL_LOOP(cgh) \ - if (true) { \ - int _rll; \ - cGH * const _cgh = const_cast<cGH*>(cgh); \ - assert (reflevel==-1); \ - for (int _rl=0; _rl<hh->reflevels(); ++_rl) { \ - set_reflevel (_cgh, _rl); \ - { -#define END_REFLEVEL_LOOP \ - } \ - } \ - set_reflevel (_cgh, -1); \ - _rll = 0; \ - } else - - - - // Reverse refinement level iterator - -#define BEGIN_REVERSE_REFLEVEL_LOOP(cgh) \ - if (true) { \ - int _rrll; \ - cGH * const _cgh = const_cast<cGH*>(cgh); \ - assert (reflevel==-1); \ - for (int _rl=hh->reflevels()-1; _rl>=0; --_rl) { \ - set_reflevel (_cgh, _rl); \ - { -#define END_REVERSE_REFLEVEL_LOOP \ - } \ - } \ - set_reflevel (_cgh, -1); \ - _rrll = 0; \ - } else - - - - // Multigrid level iterator - -#define BEGIN_MGLEVEL_LOOP(cgh) \ - if (true) { \ - int _mgl; \ - cGH * const _cgh = const_cast<cGH*>(cgh); \ - assert (reflevel>=0 && reflevel<hh->reflevels()); \ - assert (mglevel==-1); \ - for (int _ml=mglevels-1; _ml>=0; --_ml) { \ - set_mglevel (_cgh, _ml); \ - { -#define END_MGLEVEL_LOOP \ - } \ - } \ - set_mglevel (_cgh, -1); \ - _mgl = 0; \ - } else - - - - // Component iterator - - // Loop over all components. If grouptype is CCTK_GF, then loop - // over grid function components. If grouptype is CCTK_ARRAY (or - // CCTK_SCALAR), then loop over grid array (or grid scalar) - // components. In the latter case, component denotes the current - // grid array component, i.e. it cannot be used to index grid - // functions. -#define BEGIN_COMPONENT_LOOP(cgh, grouptype) \ - if (true) { \ - int _cl; \ - cGH * const _cgh = const_cast<cGH*>(cgh); \ - int const _grouptype = (grouptype); \ - int const _savec = (component); \ - int _minc, _maxc; \ - if (_grouptype == CCTK_GF) { \ - assert (reflevel>=0 && reflevel<hh->reflevels()); \ - assert (mglevel>=0 && mglevel<mglevels); \ - assert (component==-1); \ - _minc=0; \ - _maxc=hh->components(reflevel); \ - } else { \ - _minc=0; \ - _maxc=CCTK_nProcs(_cgh); \ - } \ - for (int _c=_minc; _c<_maxc; ++_c) { \ - if (component!=_c) set_component (_cgh, _c); \ - { -#define END_COMPONENT_LOOP \ - } \ - } \ - if (component!=_savec) set_component (_cgh, -1); \ - _cl = 0; \ - } else - - - - // Loop over all processor-local components. If grouptype is - // CCTK_GF, then loop over grid function components. If grouptype - // is CCTK_ARRAY (or CCTK_SCALAR), then loop over grid array (or - // grid scalar) components, i.e. execute the loop just once. In the - // latter case, component denotes the current grid array component, - // i.e. it cannot be used to index grid functions. -#define BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) \ - if (true) { \ - int _lcl; \ - cGH * const _cgh = const_cast<cGH*>(cgh); \ - int const _grouptype = (grouptype); \ - int const _savec = component; \ - int _minc, _maxc; \ - if (_grouptype == CCTK_GF) { \ - assert (reflevel>=0 && reflevel<hh->reflevels()); \ - assert (mglevel>=0 && mglevel<mglevels); \ - assert (hh->local_components(reflevel)==1 || component==-1); \ - _minc=0; \ - _maxc=hh->components(reflevel); \ - } else { \ - _minc=component; \ - _maxc=_minc+1; \ - } \ - for (int _c=_minc; _c<_maxc; ++_c) { \ - if (_grouptype!=CCTK_GF || hh->is_local(reflevel,_c)) { \ - if (component!=_c) set_component (_cgh, _c); \ - { -#define END_LOCAL_COMPONENT_LOOP \ - } \ - } \ - } \ - if (component!=_savec) set_component (_cgh, _savec); \ - _lcl = 0; \ - } else - -} // namespace Carpet +// Other declarations +#include "defines.hh" +#include "functions.hh" +#include "modes.hh" +#include "variables.hh" #endif // !defined(CARPET_PUBLIC_HH) diff --git a/Carpet/Carpet/src/make.code.defn b/Carpet/Carpet/src/make.code.defn index 893ef470d..548da3678 100644 --- a/Carpet/Carpet/src/make.code.defn +++ b/Carpet/Carpet/src/make.code.defn @@ -1,5 +1,5 @@ # Main make.code.defn file for thorn Carpet -*-Makefile-*- -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/make.code.defn,v 1.5 2002/03/23 20:20:55 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/make.code.defn,v 1.6 2004/01/25 14:57:28 schnetter Exp $ # Source files in this directory SRCS = CallFunction.cc \ @@ -17,6 +17,7 @@ SRCS = CallFunction.cc \ Shutdown.cc \ Storage.cc \ helpers.cc \ + modes.cc \ variables.cc # Subdirectories containing source files diff --git a/Carpet/CarpetIOASCII/README b/Carpet/CarpetIOASCII/README index d91cac33b..b47824330 100644 --- a/Carpet/CarpetIOASCII/README +++ b/Carpet/CarpetIOASCII/README @@ -1,7 +1,8 @@ Cactus Code Thorn CarpetIOASCII Authors : Erik Schnetter <schnetter@uni-tuebingen.de> -CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/README,v 1.2 2002/03/26 15:34:20 schnetter Exp $ +CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/README,v 1.3 2004/01/25 14:57:28 schnetter Exp $ -------------------------------------------------------------------------- Purpose of the thorn: +This thorn provides ASCII output for Carpet. diff --git a/Carpet/CarpetInterp/README b/Carpet/CarpetInterp/README index b23a9c695..845afc516 100644 --- a/Carpet/CarpetInterp/README +++ b/Carpet/CarpetInterp/README @@ -1,4 +1,4 @@ -CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/README,v 1.1 2003/04/29 14:02:25 schnetter Exp $ +CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/README,v 1.2 2004/01/25 14:57:29 schnetter Exp $ Cactus Code Thorn CarpetInterp Thorn Author(s) : Erik Schnetter <schnetter@uni-tuebingen.de> @@ -7,3 +7,4 @@ Thorn Maintainer(s) : Erik Schnetter <schnetter@uni-tuebingen.de> Purpose of the thorn: +This thorn provides a parallel interpolator for Carpet. diff --git a/Carpet/CarpetLib/README b/Carpet/CarpetLib/README index 6ee7fb750..050a838b0 100644 --- a/Carpet/CarpetLib/README +++ b/Carpet/CarpetLib/README @@ -1,7 +1,8 @@ Cactus Code Thorn CarpetLib Authors : Erik Schnetter <schnetter@uni-tuebingen.de> -CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/README,v 1.2 2002/03/26 15:34:21 schnetter Exp $ +CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/README,v 1.3 2004/01/25 14:57:29 schnetter Exp $ -------------------------------------------------------------------------- Purpose of the thorn: +This thorn contains the backend library that provides mesh refinement. diff --git a/Carpet/CarpetRegrid/schedule.ccl b/Carpet/CarpetRegrid/schedule.ccl index bdf2b7238..a0479a01d 100644 --- a/Carpet/CarpetRegrid/schedule.ccl +++ b/Carpet/CarpetRegrid/schedule.ccl @@ -1,13 +1,8 @@ # Schedule definitions for thorn CarpetRegrid -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/schedule.ccl,v 1.5 2003/09/11 16:04:15 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/schedule.ccl,v 1.6 2004/01/25 14:57:30 schnetter Exp $ schedule CarpetRegridParamcheck at PARAMCHECK { LANG: C -} "Check Parameters" - -schedule CarpetRegridStartup at STARTUP -{ - LANG: C OPTIONS: global -} "Register the regridding routine" +} "Check Parameters" diff --git a/Carpet/CarpetSlab/README b/Carpet/CarpetSlab/README index 784d30425..ea6ac8893 100644 --- a/Carpet/CarpetSlab/README +++ b/Carpet/CarpetSlab/README @@ -1,7 +1,8 @@ Cactus Code Thorn CarpetSlab Authors : Erik Schnetter <schnetter@uni-tuebingen.de> -CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/README,v 1.1 2001/03/01 13:40:10 eschnett Exp $ +CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/README,v 1.2 2004/01/25 14:57:30 schnetter Exp $ -------------------------------------------------------------------------- Purpose of the thorn: +This thorn provides hyperslabbing for Carpet. diff --git a/Carpet/CarpetSlab/src/slab.h b/Carpet/CarpetSlab/src/slab.h index 221ecf009..6e6bc3029 100644 --- a/Carpet/CarpetSlab/src/slab.h +++ b/Carpet/CarpetSlab/src/slab.h @@ -1,4 +1,4 @@ -/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.h,v 1.3 2004/01/22 13:31:07 tradke Exp $ */ +/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.h,v 1.4 2004/01/25 14:57:31 schnetter Exp $ */ #ifndef CARPETSLAB_H #define CARPETSLAB_H @@ -9,16 +9,16 @@ namespace CarpetSlab { extern "C" { #endif - - CCTK_INT CarpetSlab_Get (cGH const * const cctkGH, + + CCTK_INT CarpetSlab_Get (CCTK_POINTER_TO_CONST const cctkGH, CCTK_INT const mapping_handle, CCTK_INT const proc, CCTK_INT const vindex, CCTK_INT const timelevel, CCTK_INT const hdatatype, void * const hdata); - - CCTK_INT CarpetSlab_GetList (cGH const * const cctkGH, + + CCTK_INT CarpetSlab_GetList (CCTK_POINTER_TO_CONST const cctkGH, CCTK_INT const mapping_handle, CCTK_INT const num_arrays, CCTK_INT const * const procs, @@ -27,8 +27,8 @@ namespace CarpetSlab { CCTK_INT const * const hdatatypes, void * const * const hdata, CCTK_INT * const retvals); - - CCTK_INT CarpetSlab_LocalMappingByIndex (cGH const * const cctkGH, + + CCTK_INT CarpetSlab_LocalMappingByIndex (CCTK_POINTER_TO_CONST const cctkGH, CCTK_INT const vindex, CCTK_INT const hdim, CCTK_INT const * const direction, @@ -46,8 +46,8 @@ namespace CarpetSlab { CCTK_INT * const hsize_local, CCTK_INT * const hsize_global, CCTK_INT * const hoffset_global); - - CCTK_INT CarpetSlab_GlobalMappingByIndex (cGH const * const cctkGH, + + CCTK_INT CarpetSlab_GlobalMappingByIndex (CCTK_POINTER_TO_CONST const cctkGH, CCTK_INT const vindex, CCTK_INT const hdim, CCTK_INT const * const direction, @@ -63,13 +63,13 @@ namespace CarpetSlab { void const * const from, void * const to), CCTK_INT * const hsize); - + CCTK_INT CarpetSlab_FreeMapping (CCTK_INT const mapping_handle); /* Old interface -- don't use */ - int Hyperslab_GetHyperslab (cGH* const GH, + int Hyperslab_GetHyperslab (const cGH* const GH, const int target_proc, const int vindex, const int vtimelvl, diff --git a/Carpet/CarpetTest/src/carpettest_check_arguments.F77 b/Carpet/CarpetTest/src/carpettest_check_arguments.F77 index e4b0379fa..297e38977 100644 --- a/Carpet/CarpetTest/src/carpettest_check_arguments.F77 +++ b/Carpet/CarpetTest/src/carpettest_check_arguments.F77 @@ -1,13 +1,15 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetTest/src/carpettest_check_arguments.F77,v 1.4 2001/11/02 10:59:02 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetTest/src/carpettest_check_arguments.F77,v 1.5 2004/01/25 14:57:31 schnetter Exp $ #include "cctk.h" #include "cctk_Arguments.h" +#include "cctk_Functions.h" #include "cctk_Parameters.h" subroutine carpettest_check_arguments (CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS DECLARE_CCTK_PARAMETERS integer i,j,k print *, "Xgfg ", Xgfg0, Xgfg1, Xgfg2 @@ -16,16 +18,16 @@ c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetTest/src/carpettest_c print *, "Xarrg1 ", Xarrg10 print *, "Xscg" print * - do k=1,XGfg2 - do j=1,XGfg1 - do i=1,XGfg0 + do k=1,Xgfg2 + do j=1,Xgfg1 + do i=1,Xgfg0 gf(i,j,k) = i*10000 + j*100 + k end do end do end do - do k=1,XArrg32 - do j=1,XArrg31 - do i=1,XArrg30 + do k=1,Xarrg32 + do j=1,Xarrg31 + do i=1,Xarrg30 arr3(i,j,k) = i*10000 + j*100 + k end do end do diff --git a/CarpetAttic/CarpetIOFlexIO/README b/CarpetAttic/CarpetIOFlexIO/README index 63057198a..77a46db46 100644 --- a/CarpetAttic/CarpetIOFlexIO/README +++ b/CarpetAttic/CarpetIOFlexIO/README @@ -1,8 +1,8 @@ Cactus Code Thorn CarpetIOFlexIO Authors : Erik Schnetter <schnetter@uni-tuebingen.de> -CVS info : $Header: /home/eschnett/C/carpet/Carpet/CarpetAttic/CarpetIOFlexIO/README,v 1.3 2002/12/12 12:57:07 schnetter Exp $ +CVS info : $Header: /home/eschnett/C/carpet/Carpet/CarpetAttic/CarpetIOFlexIO/README,v 1.4 2004/01/25 14:57:29 schnetter Exp $ -------------------------------------------------------------------------- Purpose of the thorn: -FlexIO based file I/O +This thorn provides FlexIO based file I/O for Carpet. diff --git a/CarpetDev/CarpetCG/src/CG.cc b/CarpetDev/CarpetCG/src/CG.cc index 37d45e8d6..e27e465be 100644 --- a/CarpetDev/CarpetCG/src/CG.cc +++ b/CarpetDev/CarpetCG/src/CG.cc @@ -1,4 +1,4 @@ -/* $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/src/CG.cc,v 1.3 2003/11/20 08:34:03 hawke Exp $ */ +/* $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/src/CG.cc,v 1.4 2004/01/25 14:57:28 schnetter Exp $ */ #include <cassert> #include <cmath> @@ -680,345 +680,280 @@ namespace CarpetCG { assert (minerror >= 0); assert (sminerror >= 0); assert (stepsize > 0); - - // Fill common block - common::var = var; - common::res = res; - common::nvars = nvars; - common::options_table = options_table; - common::calcres = calcres; - common::applybnds = applybounds; - common::userdata = userdata; - - // Check state - // (want global mode, or at least level mode) - assert (reflevel==-1 || component==-1); - const int saved_reflevel = reflevel; - const int saved_mglevel = mglevel; - if (reflevel!=-1) { - set_mglevel ((cGH *)cctkGH, -1); - set_reflevel ((cGH *)cctkGH, -1); - } - - - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, - "Solving through nonlinear conjugate gradient iterations"); - } - - - - /* - * Literature: - * - * Jonathan Richard Shewchuk, "An Introduction to the Conjugate - * Gradient Method Without the Agonizing Pain", Technical Report - * CMU-CS-94-125, Aug. 1994 - * - * Available from http://www-2.cs.cmu.edu/~jrs/jrspapers.html - */ + + // Switch to global mode + BEGIN_GLOBAL_MODE(cctkGH) { - /* - * Notation: - * - * here there meaning - * - * iter i current cg iteration - * siter j current secant iteration - * iter2 k current cg iteration since last restart - * maxiters2 n PARA cg iterations before restart - * alpha secant step size - * beta cg step size - * delta cg distance - * epsilon PARA secant error tolerance (straight epsilon) - * minerror vareps PARA cg error tolerance (curved epsilon) - * eta secant distance - * stepsize sigma PARA secand method step parameter - * - * dir d cg direction - * -calcres f' PARA nonlinear operator - * res r cg residual - * prs s preconditioned cg residual - * var x INIT unknown variable - * - * -wgt M preconditioning matrix (not in this version) - */ + // Fill common block + common::var = var; + common::res = res; + common::nvars = nvars; + common::options_table = options_table; + common::calcres = calcres; + common::applybnds = applybounds; + common::userdata = userdata; - /* - * Algorithm: - * (Preconditioned nonlinear conjugate gradients with secant and - * Polak-Ribière) - * - * 01. i <= 0 - * 02. k <= 0 - * 03. r <= - f'(x) - * 04. calculate M \approx f''(x) - * 05. s <= M^-1 r - * 06. d <= s - * 07. delta_new <= r^T d - * 08. delta_0 <= delta_new - * 09. WHILE i < i_max AND delta_new > vareps^2 delta_0 DO - * 10. j <= 0 - * 11. delta_d <= d^T d - * 12. alpha <= - sigma_0 - * 13. eta_prev <= [f'(x + sigma_0 d)]^T d - * 14. DO - * 15. eta <= [f'(x)]^T d - * 16. alpha <= alpha (eta) / (eta_prev - eta) - * 17. x <= x + alpha d - * 18. eta_prev <= eta - * 19. j <= j + 1 - * 20. WHILE j < j_max AND alpha^2 delta_d > epsilon^2 - * 21. r <= - f'(x) - * 22. delta_old <= delta_new - * 23. delta_mid <= r^T s - * 24. calculate M \approx f''(x) - * 25. s <= M^-1 r - * 26. delta_new <= r^T s - * 27. beta <= (delta_new - delta_mid) / (delta_old) - * 28. k <= k + 1 - * 29. IF k = n OR beta <= 0 THEN - * 30. d <= s - * 31. k <= 0 - * 32. ELSE - * 33. d <= s + beta d - * 34. i <= i + 1 - */ + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Solving through nonlinear conjugate gradient iterations"); + } + + + + /* + * Literature: + * + * Jonathan Richard Shewchuk, "An Introduction to the Conjugate + * Gradient Method Without the Agonizing Pain", Technical Report + * CMU-CS-94-125, Aug. 1994 + * + * Available from http://www-2.cs.cmu.edu/~jrs/jrspapers.html + */ + + /* + * Notation: + * + * here there meaning + * + * iter i current cg iteration + * siter j current secant iteration + * iter2 k current cg iteration since last restart + * maxiters2 n PARA cg iterations before restart + * alpha secant step size + * beta cg step size + * delta cg distance + * epsilon PARA secant error tolerance (straight epsilon) + * minerror vareps PARA cg error tolerance (curved epsilon) + * eta secant distance + * stepsize sigma PARA secand method step parameter + * + * dir d cg direction + * -calcres f' PARA nonlinear operator + * res r cg residual + * prs s preconditioned cg residual + * var x INIT unknown variable + * + * -wgt M preconditioning matrix (not in this version) + */ + + /* + * Algorithm: + * (Preconditioned nonlinear conjugate gradients with secant and + * Polak-Ribière) + * + * 01. i <= 0 + * 02. k <= 0 + * 03. r <= - f'(x) + * 04. calculate M \approx f''(x) + * 05. s <= M^-1 r + * 06. d <= s + * 07. delta_new <= r^T d + * 08. delta_0 <= delta_new + * 09. WHILE i < i_max AND delta_new > vareps^2 delta_0 DO + * 10. j <= 0 + * 11. delta_d <= d^T d + * 12. alpha <= - sigma_0 + * 13. eta_prev <= [f'(x + sigma_0 d)]^T d + * 14. DO + * 15. eta <= [f'(x)]^T d + * 16. alpha <= alpha (eta) / (eta_prev - eta) + * 17. x <= x + alpha d + * 18. eta_prev <= eta + * 19. j <= j + 1 + * 20. WHILE j < j_max AND alpha^2 delta_d > epsilon^2 + * 21. r <= - f'(x) + * 22. delta_old <= delta_new + * 23. delta_mid <= r^T s + * 24. calculate M \approx f''(x) + * 25. s <= M^-1 r + * 26. delta_new <= r^T s + * 27. beta <= (delta_new - delta_mid) / (delta_old) + * 28. k <= k + 1 + * 29. IF k = n OR beta <= 0 THEN + * 30. d <= s + * 31. k <= 0 + * 32. ELSE + * 33. d <= s + beta d + * 34. i <= i + 1 + */ - // Pointers to grid variables and temporary storage - int varptrs[nvars]; - int resptrs[nvars]; - // CCTK_REAL * restrict * restrict prsptrs; - // CCTK_REAL * restrict * restrict dirptrs; + // Pointers to grid variables and temporary storage + vector<int> varptrs(nvars); + vector<int> resptrs(nvars); + // CCTK_REAL * restrict * restrict prsptrs; + // CCTK_REAL * restrict * restrict dirptrs; - // Get storage pointers - for (n=0; n<nvars; ++n) { - varptrs[n] = var[n]; - resptrs[n] = res[n]; - } + // Get storage pointers + for (n=0; n<nvars; ++n) { + varptrs[n] = var[n]; + resptrs[n] = res[n]; + } - // Allocate temporary memory - // prsptrs = malloc (nvars * sizeof *prsptrs); - // dirptrs = malloc (nvars * sizeof *prsptrs); - assert (cg_maxsolvevars >= nvars); - int prsptrs[cg_maxsolvevars]; - int dirptrs[cg_maxsolvevars]; - assert (prsptrs); - assert (dirptrs); - for (n=0; n<nvars; ++n) { - char varname[1000]; - sprintf(varname,"CarpetCG::prsvars[%d]",n); - prsptrs[n] = CCTK_VarIndex(varname); - assert ((prsptrs[n] >= 0)&&(prsptrs[n] < CCTK_NumVars())); - sprintf(varname,"CarpetCG::dirvars[%d]",n); - dirptrs[n] = CCTK_VarIndex(varname); - assert ((dirptrs[n] >= 0)&&(dirptrs[n] < CCTK_NumVars())); - } + // Allocate temporary memory + // prsptrs = malloc (nvars * sizeof *prsptrs); + // dirptrs = malloc (nvars * sizeof *prsptrs); + assert (cg_maxsolvevars >= nvars); + vector<int> prsptrs(cg_maxsolvevars); + vector<int> dirptrs(cg_maxsolvevars); + for (n=0; n<nvars; ++n) { + char varname[1000]; + sprintf(varname,"CarpetCG::prsvars[%d]",n); + prsptrs[n] = CCTK_VarIndex(varname); + assert ((prsptrs[n] >= 0)&&(prsptrs[n] < CCTK_NumVars())); + sprintf(varname,"CarpetCG::dirvars[%d]",n); + dirptrs[n] = CCTK_VarIndex(varname); + assert ((dirptrs[n] >= 0)&&(dirptrs[n] < CCTK_NumVars())); + } - nexttime = time(0); + nexttime = time(0); - /* 01. i <= 0 */ - iter = 0; + /* 01. i <= 0 */ + iter = 0; - /* 02. k <= 0 */ - iter2 = 0; + /* 02. k <= 0 */ + iter2 = 0; - /* 03. r <= - f'(x) */ - /* 04. calculate M \approx f''(x) */ - common::ierr = 0; - CallLocalFunction((cGH *)cctkGH, common::call_calcres); - ierr = common::ierr; - assert (! ierr); - - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - } - CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); - - /* 05. s <= M^-1 r */ - // No preconditioning - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - common::toindex[n] = prsptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } - - // cout << "dirptrs " << dirptrs[0] << endl; - - CallLocalFunction((cGH *)cctkGH, common::call_apply_preconditioner); - - // cout << "dirptrs " << dirptrs[0] << endl; - - /* 06. d <= s */ - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = prsptrs[n]; - common::toindex[n] = dirptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } - - // cout << "dirptrs " << dirptrs[0] << endl; - - CallLocalFunction((cGH *)cctkGH, common::call_copy); - - // cout << "dirptrs " << dirptrs[0] << endl; - - /* 07. delta_new <= r^T d */ - // output (cctkGH, var, res, wgt, nvars, iter); - common::realoutput_count = 0; - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = prsptrs[n]; - common::toindex[n] = resptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } - - // cout << "dirptrs " << dirptrs[0] << endl; - - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - delta_new = common::realoutput; - gsize = common::realoutput_count; - // cout << "delta_new " << delta_new << " gsize " << gsize << endl; - - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - common::toindex[n] = resptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } - - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - epsilon = common::realoutput; -// cout << "epsilon " << epsilon << endl; - - /* 08. delta_0 <= delta_new */ - delta_0 = delta_new; - - /* 09. WHILE i < i_max AND delta_new > vareps^2 delta_0 DO */ - while (iter < maxiters && epsilon / (nvars * gsize) > pow(minerror,2)) { - - if (verbose || veryverbose) { - currenttime = time(0); - if (veryverbose || (iter % outevery == 0 && currenttime >= nexttime)) { - CCTK_VInfo (CCTK_THORNSTRING, - "Iteration %d (%d since restart): residual is %g (%g,%d,%g)", - iter, iter2, (double)sqrt(epsilon / (nvars * gsize)),epsilon,nvars,gsize); - if (outeveryseconds > 0) { - while (nexttime <= currenttime) nexttime += outeveryseconds; - } - } - } - - /* 10. j <= 0 */ - siter = 0; + /* 03. r <= - f'(x) */ + /* 04. calculate M \approx f''(x) */ + common::ierr = 0; + CallLocalFunction((cGH *)cctkGH, common::call_calcres); + ierr = common::ierr; + assert (! ierr); - /* 11. delta_d <= d^T d */ for (int n = 0; n < nvars; n++) { - common::fromindex[n] = dirptrs[n]; - common::toindex[n] = dirptrs[n]; + common::fromindex[n] = resptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; - common::toindex[n] = -1; } - - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - delta_d = common::realoutput; - -// cout << "delta_d " << delta_d << endl; - - /* 12. alpha <= - sigma_0 */ - alpha = - stepsize; - sum_alpha = alpha; - do_abort = 0; + CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); - /* 13. eta_prev <= [f'(x + sigma_0 d)]^T d */ + /* 05. s <= M^-1 r */ + // No preconditioning for (int n = 0; n < nvars; n++) { - common::fromindex[n] = dirptrs[n]; - common::toindex[n] = varptrs[n]; + common::fromindex[n] = resptrs[n]; + common::toindex[n] = prsptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; common::toindex[n] = -1; } - common::realconstant = -alpha; -// cout << "Add scaled, const " << -alpha << endl; - - CallLocalFunction((cGH *)cctkGH, common::call_add_scaled); - common::ierr = 0; - CallLocalFunction((cGH *)cctkGH, common::call_applybnds); - assert (! ierr); - ierr = 0; - CallLocalFunction((cGH *)cctkGH, common::call_calcres); - assert (! ierr); - + // cout << "dirptrs " << dirptrs[0] << endl; + + CallLocalFunction((cGH *)cctkGH, common::call_apply_preconditioner); + + // cout << "dirptrs " << dirptrs[0] << endl; + + /* 06. d <= s */ for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; + common::fromindex[n] = prsptrs[n]; + common::toindex[n] = dirptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; + common::toindex[n] = -1; } - CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); + // cout << "dirptrs " << dirptrs[0] << endl; + + CallLocalFunction((cGH *)cctkGH, common::call_copy); + + // cout << "dirptrs " << dirptrs[0] << endl; + + /* 07. delta_new <= r^T d */ + // output (cctkGH, var, res, wgt, nvars, iter); + common::realoutput_count = 0; for (int n = 0; n < nvars; n++) { - common::fromindex[n] = dirptrs[n]; + common::fromindex[n] = prsptrs[n]; common::toindex[n] = resptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; common::toindex[n] = -1; } + + // cout << "dirptrs " << dirptrs[0] << endl; + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - eta = - common::realoutput; - -// cout << "eta " << eta << endl; - + delta_new = common::realoutput; + gsize = common::realoutput_count; + // cout << "delta_new " << delta_new << " gsize " << gsize << endl; + for (int n = 0; n < nvars; n++) { - common::fromindex[n] = dirptrs[n]; - common::toindex[n] = varptrs[n]; + common::fromindex[n] = resptrs[n]; + common::toindex[n] = resptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; common::toindex[n] = -1; } - common::realconstant = alpha; -// cout << "Add scaled, const " << alpha << endl; - - CallLocalFunction((cGH *)cctkGH, common::call_add_scaled); - ierr = 0; - CallLocalFunction((cGH *)cctkGH, common::call_applybnds); - assert (! ierr); - eta_prev = eta; + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + epsilon = common::realoutput; + // cout << "epsilon " << epsilon << endl; - /* 14. DO */ - do { + /* 08. delta_0 <= delta_new */ + delta_0 = delta_new; + + /* 09. WHILE i < i_max AND delta_new > vareps^2 delta_0 DO */ + while (iter < maxiters && epsilon / (nvars * gsize) > pow(minerror,2)) { - if (veryverbose) { - CCTK_VInfo - (CCTK_THORNSTRING, - " Secant iteration %d: step size is %g, orthogonality is %g", - siter, - (double)alpha, (double)sqrt(fabs(eta_prev) / (nvars * gsize))); + if (verbose || veryverbose) { + currenttime = time(0); + if (veryverbose || (iter % outevery == 0 && currenttime >= nexttime)) { + CCTK_VInfo (CCTK_THORNSTRING, + "Iteration %d (%d since restart): residual is %g (%g,%d,%g)", + iter, iter2, (double)sqrt(epsilon / (nvars * gsize)),epsilon,nvars,gsize); + if (outeveryseconds > 0) { + while (nexttime <= currenttime) nexttime += outeveryseconds; + } + } } - - /* 15. eta <= [f'(x)]^T d */ + + /* 10. j <= 0 */ + siter = 0; + + /* 11. delta_d <= d^T d */ + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = dirptrs[n]; + common::toindex[n] = dirptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } + + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + delta_d = common::realoutput; + + // cout << "delta_d " << delta_d << endl; + + /* 12. alpha <= - sigma_0 */ + alpha = - stepsize; + sum_alpha = alpha; + do_abort = 0; + + /* 13. eta_prev <= [f'(x + sigma_0 d)]^T d */ + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = dirptrs[n]; + common::toindex[n] = varptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } + common::realconstant = -alpha; + // cout << "Add scaled, const " << -alpha << endl; + + CallLocalFunction((cGH *)cctkGH, common::call_add_scaled); + common::ierr = 0; + CallLocalFunction((cGH *)cctkGH, common::call_applybnds); + assert (! ierr); + ierr = 0; CallLocalFunction((cGH *)cctkGH, common::call_calcres); assert (! ierr); @@ -1030,7 +965,7 @@ namespace CarpetCG { common::fromindex[n] = -1; } CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); - + for (int n = 0; n < nvars; n++) { common::fromindex[n] = dirptrs[n]; common::toindex[n] = resptrs[n]; @@ -1039,34 +974,10 @@ namespace CarpetCG { common::fromindex[n] = -1; common::toindex[n] = -1; } - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); eta = - common::realoutput; -// cout << "eta " << eta << endl; - - /* 16. alpha <= alpha (eta) / (eta_prev - eta) */ - alpha_old = alpha; - alpha *= eta / (eta_prev - eta); - sum_alpha += alpha; - if (veryverbose) { - CCTK_VInfo - (CCTK_THORNSTRING, - " Changing step size, iteration %d: was %g, now %g (%g, %g)", - siter, alpha_old, alpha, eta, eta_prev); - } - assert (sum_alpha > - maxstepsize); - if (sum_alpha > maxstepsize) { - alpha = maxstepsize - alpha_old; - do_abort = 1; - if (verbose) { - CCTK_VInfo - (CCTK_THORNSTRING, - " Secant iteration %d: limiting total step size", siter); - } - } - - /* 17. x <= x + alpha d */ + // cout << "eta " << eta << endl; for (int n = 0; n < nvars; n++) { common::fromindex[n] = dirptrs[n]; @@ -1077,174 +988,249 @@ namespace CarpetCG { common::toindex[n] = -1; } common::realconstant = alpha; -// cout << "Add scaled, const " << alpha << endl; + // cout << "Add scaled, const " << alpha << endl; CallLocalFunction((cGH *)cctkGH, common::call_add_scaled); ierr = 0; CallLocalFunction((cGH *)cctkGH, common::call_applybnds); assert (! ierr); - - /* 18. eta_prev <= eta */ + eta_prev = eta; - - /* 19. j <= j + 1 */ - ++ siter; - - /* 20. WHILE j < j_max AND alpha^2 delta_d > epsilon^2 */ - } while (siter < smaxiters - && pow(alpha,2) * delta_d > pow(sminerror,2) - && !do_abort); - if (veryverbose) { - CCTK_VInfo - (CCTK_THORNSTRING, - " Secant iteration %d: step size is %g, orthogonality is %g", - siter, - (double)alpha, (double)sqrt(fabs(eta_prev) / (nvars * gsize))); - } + /* 14. DO */ + do { - /* 21. r <= - f'(x) */ - /* 24. calculate M \approx f''(x) */ - ierr = 0; - CallLocalFunction((cGH *)cctkGH, common::call_calcres); - assert (! ierr); + if (veryverbose) { + CCTK_VInfo + (CCTK_THORNSTRING, + " Secant iteration %d: step size is %g, orthogonality is %g", + siter, + (double)alpha, (double)sqrt(fabs(eta_prev) / (nvars * gsize))); + } + + /* 15. eta <= [f'(x)]^T d */ + ierr = 0; + CallLocalFunction((cGH *)cctkGH, common::call_calcres); + assert (! ierr); - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - } - CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); - - /* 23. delta_mid <= r^T s */ - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = prsptrs[n]; - common::toindex[n] = resptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = resptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + } + CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - delta_mid = common::realoutput; -// cout << "delta_mid " << delta_mid << endl; - - /* 25. s <= M^-1 r */ - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - common::toindex[n] = prsptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } - - CallLocalFunction((cGH *)cctkGH, common::call_apply_preconditioner); + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = dirptrs[n]; + common::toindex[n] = resptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } - /* 22. delta_old <= delta_new */ - delta_old = delta_new; - - /* 26. delta_new <= r^T s */ -// output (cctkGH, var, res, wgt, nvars, iter+1); - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = prsptrs[n]; - common::toindex[n] = resptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + eta = - common::realoutput; - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - delta_new = common::realoutput; -// cout << "delta_new " << delta_new << endl; + // cout << "eta " << eta << endl; + + /* 16. alpha <= alpha (eta) / (eta_prev - eta) */ + alpha_old = alpha; + alpha *= eta / (eta_prev - eta); + sum_alpha += alpha; + if (veryverbose) { + CCTK_VInfo + (CCTK_THORNSTRING, + " Changing step size, iteration %d: was %g, now %g (%g, %g)", + siter, alpha_old, alpha, eta, eta_prev); + } + assert (sum_alpha > - maxstepsize); + if (sum_alpha > maxstepsize) { + alpha = maxstepsize - alpha_old; + do_abort = 1; + if (verbose) { + CCTK_VInfo + (CCTK_THORNSTRING, + " Secant iteration %d: limiting total step size", siter); + } + } + + /* 17. x <= x + alpha d */ - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - common::toindex[n] = resptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = dirptrs[n]; + common::toindex[n] = varptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } + common::realconstant = alpha; + // cout << "Add scaled, const " << alpha << endl; - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - epsilon = common::realoutput; -// cout << "epsilon " << epsilon << endl; + CallLocalFunction((cGH *)cctkGH, common::call_add_scaled); + ierr = 0; + CallLocalFunction((cGH *)cctkGH, common::call_applybnds); + assert (! ierr); + + /* 18. eta_prev <= eta */ + eta_prev = eta; + + /* 19. j <= j + 1 */ + ++ siter; + + /* 20. WHILE j < j_max AND alpha^2 delta_d > epsilon^2 */ + } while (siter < smaxiters + && pow(alpha,2) * delta_d > pow(sminerror,2) + && !do_abort); - /* 27. beta <= (delta_new - delta_mid) / (delta_old) */ - beta = (delta_new - delta_mid) / delta_old; + if (veryverbose) { + CCTK_VInfo + (CCTK_THORNSTRING, + " Secant iteration %d: step size is %g, orthogonality is %g", + siter, + (double)alpha, (double)sqrt(fabs(eta_prev) / (nvars * gsize))); + } - /* 28. k <= k + 1 */ - ++ iter2; + /* 21. r <= - f'(x) */ + /* 24. calculate M \approx f''(x) */ + ierr = 0; + CallLocalFunction((cGH *)cctkGH, common::call_calcres); + assert (! ierr); + + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = resptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + } + CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); - /* 29. IF k = n OR beta <= 0 THEN */ - if (iter2 >= maxiters2 || beta <= 0 || do_abort) { - /* Restart */ - - /* 30. d <= s */ + /* 23. delta_mid <= r^T s */ for (int n = 0; n < nvars; n++) { common::fromindex[n] = prsptrs[n]; - common::toindex[n] = dirptrs[n]; + common::toindex[n] = resptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; common::toindex[n] = -1; } - CallLocalFunction((cGH *)cctkGH, common::call_copy); - - /* 31. k <= 0 */ - iter2 = 0; + + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + delta_mid = common::realoutput; + // cout << "delta_mid " << delta_mid << endl; - /* 32. ELSE */ - } else { - /* Continue with CG iterations */ + /* 25. s <= M^-1 r */ + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = resptrs[n]; + common::toindex[n] = prsptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } - /* 33. d <= s + beta d */ + CallLocalFunction((cGH *)cctkGH, common::call_apply_preconditioner); + + /* 22. delta_old <= delta_new */ + delta_old = delta_new; + + /* 26. delta_new <= r^T s */ + // output (cctkGH, var, res, wgt, nvars, iter+1); for (int n = 0; n < nvars; n++) { common::fromindex[n] = prsptrs[n]; - common::toindex[n] = dirptrs[n]; + common::toindex[n] = resptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } + + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + delta_new = common::realoutput; + // cout << "delta_new " << delta_new << endl; + + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = resptrs[n]; + common::toindex[n] = resptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; common::toindex[n] = -1; } - common::realconstant = beta; - // cout << "beta " << beta << endl; + + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + epsilon = common::realoutput; + // cout << "epsilon " << epsilon << endl; + + /* 27. beta <= (delta_new - delta_mid) / (delta_old) */ + beta = (delta_new - delta_mid) / delta_old; + + /* 28. k <= k + 1 */ + ++ iter2; + + /* 29. IF k = n OR beta <= 0 THEN */ + if (iter2 >= maxiters2 || beta <= 0 || do_abort) { + /* Restart */ + + /* 30. d <= s */ + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = prsptrs[n]; + common::toindex[n] = dirptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } + CallLocalFunction((cGH *)cctkGH, common::call_copy); + + /* 31. k <= 0 */ + iter2 = 0; + + /* 32. ELSE */ + } else { + /* Continue with CG iterations */ + + /* 33. d <= s + beta d */ + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = prsptrs[n]; + common::toindex[n] = dirptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } + common::realconstant = beta; + // cout << "beta " << beta << endl; - CallLocalFunction((cGH *)cctkGH, common::call_scale_and_add); + CallLocalFunction((cGH *)cctkGH, common::call_scale_and_add); - } /* if iter2 */ + } /* if iter2 */ /* 34. i <= i + 1 */ - ++ iter; + ++ iter; - } /* for iter */ + } /* for iter */ - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, - "Iteration %d (%d since restart): residual is %g", - iter, iter2, (double)sqrt(epsilon / (nvars * gsize))); - } + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Iteration %d (%d since restart): residual is %g", + iter, iter2, (double)sqrt(epsilon / (nvars * gsize))); + } /* Free temporary memory */ -// for (n=0; n<nvars; ++n) { -// free (prsptrs[n]); -// free (dirptrs[n]); -// } - -// free (varptrs); -// free (resptrs); -// free (prsptrs); -// free (dirptrs); - - // Restore state - if (reflevel!=saved_reflevel) { - set_reflevel ((cGH *)cctkGH, saved_reflevel); - set_mglevel ((cGH *)cctkGH, saved_mglevel); - } + // for (n=0; n<nvars; ++n) { + // free (prsptrs[n]); + // free (dirptrs[n]); + // } + + // free (varptrs); + // free (resptrs); + // free (prsptrs); + // free (dirptrs); + + } END_GLOBAL_MODE; return 0; } diff --git a/CarpetDev/CarpetJacobi/param.ccl b/CarpetDev/CarpetJacobi/param.ccl index ce69dbc91..bdeecf41b 100644 --- a/CarpetDev/CarpetJacobi/param.ccl +++ b/CarpetDev/CarpetJacobi/param.ccl @@ -1,5 +1,5 @@ # Parameter definitions for thorn CarpetJacobi -# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/param.ccl,v 1.1 2003/09/02 14:35:57 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/param.ccl,v 1.2 2004/01/25 14:57:29 schnetter Exp $ BOOLEAN verbose "Produce log output while running" STEERABLE=always { @@ -18,3 +18,15 @@ CCTK_INT outeveryseconds "Produce output every n seconds" STEERABLE=always { 0:* :: "" } 10 + + + +CCTK_INT solve_minlevel "First level on which to solve" STEERABLE=always +{ + 0:* :: "" +} 0 + +CCTK_INT reflevelpower "Power of the refinement factor" STEERABLE=always +{ + 0:* :: "" +} 2 diff --git a/CarpetDev/CarpetJacobi/src/Jacobi.cc b/CarpetDev/CarpetJacobi/src/Jacobi.cc index ad673b342..c326b66cd 100644 --- a/CarpetDev/CarpetJacobi/src/Jacobi.cc +++ b/CarpetDev/CarpetJacobi/src/Jacobi.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/Jacobi.cc,v 1.2 2003/11/19 14:05:14 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/Jacobi.cc,v 1.3 2004/01/25 14:57:29 schnetter Exp $ #include <cassert> #include <cmath> @@ -14,14 +14,18 @@ #include "util_ErrorCodes.h" #include "util_Table.h" -#include "carpet.hh" #include "TATelliptic.h" +#include "carpet.hh" + +#include "th.hh" + namespace Carpet { // TODO: fix this - void Restrict (const cGH* cgh); + void CycleTimeLevels (const cGH* cctkGH); + void Restrict (const cGH* cctkGH); }; @@ -48,134 +52,6 @@ namespace CarpetJacobi { - namespace common { - - const int * var; - const int * res; - int nvars; - int options_table; - calcfunc calcres; - calcfunc applybnds; - void * userdata; - - CCTK_REAL factor; - CCTK_REAL * factors; - - CCTK_REAL norm_count; - CCTK_REAL norm_l2; - - int ierr; - - void call_calcres (cGH * const cctkGH) - { - if (ierr) return; - ierr = calcres (cctkGH, options_table, userdata); - } - - void norm (cGH * const cctkGH) - { - int ierr; - - cGroup groupdata; - ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata); - assert (!ierr); - cGroupDynamicData groupdyndata; - ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]), - &groupdyndata); - assert (!ierr); - const int thedim = groupdata.dim; - assert (thedim>=0 && thedim<=dim); - assert (thedim == groupdyndata.dim); - int lsh[dim], nghostzones[dim]; - for (int d=0; d<thedim; ++d) { - lsh[d] = groupdyndata.lsh[d]; - nghostzones[d] = groupdyndata.nghostzones[d]; - } - for (int d=thedim; d<dim; ++d) { - lsh[d] = 1; - nghostzones[d] = 0; - } - for (int d=0; d<dim; ++d) { - assert (lsh[d]>=0); - assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]); - } - - for (int n=0; n<nvars; ++n) { - const CCTK_REAL * restrict resptr - = (const CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, res[n]); - assert (resptr); - const CCTK_REAL fac = factors[n]; - for (int k=nghostzones[2]; k<lsh[2]-nghostzones[2]; ++k) { - for (int j=nghostzones[1]; j<lsh[1]-nghostzones[1]; ++j) { - for (int i=nghostzones[0]; i<lsh[0]-nghostzones[0]; ++i) { - const int ind = i + lsh[0] * (j + lsh[1] * k); - ++norm_count; - // TODO: scale the norm by the resolution? - norm_l2 += pow(fac * resptr[ind], 2); - } - } - } - } - } - - void smooth (cGH * const cctkGH) - { - int ierr; - - cGroup groupdata; - ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata); - assert (!ierr); - cGroupDynamicData groupdyndata; - ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]), - &groupdyndata); - assert (!ierr); - const int thedim = groupdata.dim; - assert (thedim>=0 && thedim<=dim); - assert (thedim == groupdyndata.dim); - int lsh[dim], nghostzones[dim]; - for (int d=0; d<thedim; ++d) { - lsh[d] = groupdyndata.lsh[d]; - nghostzones[d] = groupdyndata.nghostzones[d]; - } - for (int d=thedim; d<dim; ++d) { - lsh[d] = 1; - nghostzones[d] = 0; - } - for (int d=0; d<dim; ++d) { - assert (lsh[d]>=0); - assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]); - } - - for (int n=0; n<nvars; ++n) { - CCTK_REAL * restrict varptr - = (CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, var[n]); - assert (varptr); - const CCTK_REAL * restrict resptr - = (const CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, res[n]); - assert (resptr); - assert (resptr != varptr); - const CCTK_REAL fac = factor * factors[n]; - for (int k=nghostzones[2]; k<lsh[2]-nghostzones[2]; ++k) { - for (int j=nghostzones[1]; j<lsh[1]-nghostzones[1]; ++j) { - for (int i=nghostzones[0]; i<lsh[0]-nghostzones[0]; ++i) { - const int ind = i + lsh[0] * (j + lsh[1] * k); - varptr[ind] += fac * resptr[ind]; - } - } - } - } - } - - void call_applybnds (cGH * const cctkGH) - { - if (ierr) return; - ierr = applybnds (cctkGH, options_table, userdata); - } - - } // namespace common - - - void CarpetJacobi_register () { int const ierr = TATelliptic_RegisterSolver @@ -194,7 +70,6 @@ namespace CarpetJacobi { const calcfunc applybnds, void * const userdata) { - DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; // Error control @@ -212,19 +87,19 @@ namespace CarpetJacobi { assert (res); for (int n=0; n<nvars; ++n) { assert (var[n] >= 0 && var[n] < CCTK_NumVars()); - assert (CCTK_GroupTypeFromVarI(var[n]) == CCTK_GF - || CCTK_GroupTypeFromVarI(var[n]) == CCTK_ARRAY); - assert (CCTK_GroupDimFromVarI(var[n]) <= dim); + assert (CCTK_GroupTypeFromVarI(var[n]) == CCTK_GF); + assert (CCTK_GroupDimFromVarI(var[n]) == dim); assert (CCTK_VarTypeI(var[n]) == CCTK_VARIABLE_REAL); - assert (CCTK_QueryGroupStorageI(cctkGH, CCTK_GroupIndexFromVarI(var[n]))); + assert (CCTK_QueryGroupStorageI(cctkGH, + CCTK_GroupIndexFromVarI(var[n]))); } for (int n=0; n<nvars; ++n) { assert (res[n] >= 0 && res[n] < CCTK_NumVars()); - assert (CCTK_GroupTypeFromVarI(res[n]) == CCTK_GF - || CCTK_GroupTypeFromVarI(res[n]) == CCTK_ARRAY); - assert (CCTK_GroupDimFromVarI(res[n]) <= dim); + assert (CCTK_GroupTypeFromVarI(res[n]) == CCTK_GF); + assert (CCTK_GroupDimFromVarI(res[n]) == dim); assert (CCTK_VarTypeI(res[n]) == CCTK_VARIABLE_REAL); - assert (CCTK_QueryGroupStorageI(cctkGH, CCTK_GroupIndexFromVarI(res[n]))); + assert (CCTK_QueryGroupStorageI(cctkGH, + CCTK_GroupIndexFromVarI(res[n]))); } for (int n=0; n<nvars; ++n) { assert (var[n] != res[n]); @@ -236,47 +111,7 @@ namespace CarpetJacobi { } } -#if 0 - // Get domain description - cGroup groupdata; - ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata); - assert (!ierr); - cGroupDynamicData groupdyndata; - ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]), - &groupdyndata); - assert (!ierr); - const int thedim = groupdata.dim; - assert (thedim>=0 && thedim<=dim); - assert (thedim == groupdyndata.dim); - int lsh[dim], nghostzones[dim]; - for (int d=0; d<thedim; ++d) { - lsh[d] = groupdyndata.lsh[d]; - nghostzones[d] = groupdyndata.nghostzones[d]; - } - for (int d=thedim; d<dim; ++d) { - lsh[d] = 1; - nghostzones[d] = 0; - } - for (int d=0; d<dim; ++d) { - assert (lsh[d]>=0); - assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]); - } - // Check all variables - for (int n=0; n<nvars; ++n) { - ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[n]), &groupdata); - assert (!ierr); - ierr = CCTK_GroupDynamicData - (cctkGH, CCTK_GroupIndexFromVarI(var[n]), &groupdyndata); - assert (!ierr); - assert (groupdata.dim == thedim); - assert (groupdyndata.dim == thedim); - for (int d=0; d<thedim; ++d) { - assert (groupdyndata.lsh[d] == lsh[d]); - assert (groupdyndata.nghostzones[d] == nghostzones[d]); - } - } -#endif assert (options_table >= 0); @@ -310,7 +145,8 @@ namespace CarpetJacobi { assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); ierr = Util_TableGetReal (options_table, &deceleration, "deceleration"); assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); - ierr = Util_TableGetRealArray (options_table, nvars, &factors[0], "factors"); + ierr = Util_TableGetRealArray (options_table, + nvars, &factors[0], "factors"); assert (ierr==nvars || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); assert (maxiters >= 0); @@ -324,197 +160,375 @@ namespace CarpetJacobi { assert (factors[n] != 0); } - // Check state - // (want global mode, or at least level mode) - assert (reflevel==-1 || component==-1); - const int saved_reflevel = reflevel; - const int saved_mglevel = mglevel; - if (reflevel!=-1) { - set_mglevel ((cGH *)cctkGH, -1); - set_reflevel ((cGH *)cctkGH, -1); - } - // Fill common block - common::var = var; - common::res = res; - common::nvars = nvars; - common::options_table = options_table; - common::calcres = calcres; - common::applybnds = applybnds; - common::userdata = userdata; - common::factor = factor; - common::factors = &factors[0]; - // Init statistics - CCTK_REAL norm2 = HUGE_VAL; - CCTK_REAL speed = 0.0; - CCTK_REAL throttle = 1.0; - bool did_hit_min = false; - bool did_hit_max = false; - const time_t starttime = time(0); - time_t nexttime = starttime; + assert (is_global_mode() || is_level_mode()); + + // The level to solve on + const int solve_level = is_level_mode() ? reflevel : reflevels-1; - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, "Solving through adaptive Jacobi iterations"); + if (solve_level < solve_minlevel) { + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Did not solve because the level is too coarse."); + } + Util_TableSetInt (options_table, 0, "iters"); + Util_TableSetReal (options_table, -1.0, "error"); + + // Return as if the solving had not converged + return -1; } - int iter; - for (iter=0; iter<=maxiters; ++iter) { +#warning "TODO: assert that all levels are at the same time" + + // Switch to global mode + BEGIN_GLOBAL_MODE(cctkGH) { - // Calculate residual - common::ierr = 0; - CallLocalFunction ((cGH *)cctkGH, common::call_calcres); - ierr = common::ierr; - if (ierr != 0) { - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, "Residual calculation reported error; aborting solver."); + // Reset the initial data time + global_time = 0; + delta_time = 1; + * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; + * const_cast<CCTK_REAL *> (& cctkGH->cctk_delta_time) = delta_time; + BEGIN_REFLEVEL_LOOP(cctkGH) { + * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; + for (int m=0; m<maps; ++m) { + vtt[m]->set_time (reflevel, mglevel, 0); + vtt[m]->set_delta + (reflevel, mglevel, 1.0 / ipow(reflevelfact, reflevelpower)); } - Util_TableSetInt (options_table, iter, "iters"); - Util_TableDeleteKey (options_table, "error"); - goto done; - } + } END_REFLEVEL_LOOP; - // Save old norm - const CCTK_REAL old_norm2 = norm2; - // Calculate norm - common::norm_count = 0; - common::norm_l2 = 0; - CallLocalFunction ((cGH *)cctkGH, common::norm); - const int sum_handle = CCTK_ReductionArrayHandle ("sum"); - assert (sum_handle >= 0); - CCTK_REAL reduce_in[2], reduce_out[2]; - reduce_in[0] = common::norm_count; - reduce_in[1] = common::norm_l2; - ierr = CCTK_ReduceLocArrayToArray1D (cctkGH, -1, sum_handle, - reduce_in, reduce_out, 2, - CCTK_VARIABLE_REAL); - norm2 = sqrt(reduce_out[1] / reduce_out[0]); - // Calculate convergence speed - speed = old_norm2 / norm2; + // Init statistics + vector<CCTK_REAL> norm_counts (solve_level+1); + vector<CCTK_REAL> norm_l2s (solve_level+1); + CCTK_REAL norm2 = HUGE_VAL; + CCTK_REAL speed = 0.0; + CCTK_REAL throttle = 1.0; + bool did_hit_min = false; + bool did_hit_max = false; + const time_t starttime = time(0); + time_t nexttime = starttime; - // Log output if (verbose || veryverbose) { - const time_t currenttime = time(0); - if (iter == maxiters -#if HAVE_ISNAN - || isnan(norm2) -#endif - || norm2 <= minerror - || (iter % outevery == 0 && currenttime >= nexttime)) { - if (verbose || veryverbose) { - if (did_hit_min) { - CCTK_VInfo (CCTK_THORNSTRING, - "Convergence factor became too small: artificially kept at minfactor (may lead to instability)"); + CCTK_VInfo (CCTK_THORNSTRING, + "Solving through adaptive Jacobi iterations"); + } + + const int icoarsestep + = ipow (mglevelfact * maxreflevelfact, reflevelpower); + const int istep + = ipow (mglevelfact * maxreflevelfact / ipow(reffact, solve_level), + reflevelpower); + int iter = istep; + for (;; iter+=istep) { + + global_time = 1.0 * iter / ipow(maxreflevelfact, reflevelpower); + * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; +// cout << "CJ global iter " << iter << " time " << global_time << flush << endl; + + // Calculate residual, smooth, and apply boundary conditions + ierr = 0; + BEGIN_REFLEVEL_LOOP(cctkGH) { + const int do_every + = ipow (mglevelfact * (maxreflevelfact/reflevelfact), + reflevelpower); + if (reflevel <= solve_level && (iter-istep) % do_every == 0) { + + // Advance time + for (int m=0; m<maps; ++m) { + vtt[m]->advance_time (reflevel, mglevel); } - if (did_hit_max) { - CCTK_VInfo (CCTK_THORNSTRING, - "Convergence factor became too large: artificially kept at maxfactor (may lead to slower convergence)"); + * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) + = (1.0 * (iter - istep + do_every) + / ipow(maxreflevelfact, reflevelpower)); + CycleTimeLevels (cctkGH); +// cout << "CJ residual iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl; + + // Calculate residual + BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { + if (! ierr) { + ierr = calcres (cctkGH, options_table, userdata); + } + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; + + // Smooth and apply boundary conditions + CCTK_REAL levfac = 0; + for (int d=0; d<dim; ++d) { + levfac += 1 / + pow (cctkGH->cctk_delta_space[d] / cctkGH->cctk_levfac[d], 2); } - did_hit_min = false; - did_hit_max = false; + levfac = 1 / levfac; + + BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { + if (! ierr) { + + for (int n=0; n<nvars; ++n) { + CCTK_REAL * restrict varptr + = (CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, var[n]); + assert (varptr); + const CCTK_REAL * restrict resptr + = (const CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, res[n]); + assert (resptr); + assert (resptr != varptr); + const CCTK_REAL fac = factor * factors[n] * levfac; + assert (cctkGH->cctk_dim == 3); + for (int k=cctkGH->cctk_nghostzones[2]; + k<cctkGH->cctk_lsh[2]-cctkGH->cctk_nghostzones[2]; + ++k) { + for (int j=cctkGH->cctk_nghostzones[1]; + j<cctkGH->cctk_lsh[1]-cctkGH->cctk_nghostzones[1]; + ++j) { + for (int i=cctkGH->cctk_nghostzones[0]; + i<cctkGH->cctk_lsh[0]-cctkGH->cctk_nghostzones[0]; + ++i) { + const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k); + varptr[ind] += fac * resptr[ind]; + } + } + } + } // for n + + // Apply boundary conditions + ierr = applybnds (cctkGH, options_table, userdata); + + } + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; + } - if (veryverbose) { + } END_REFLEVEL_LOOP; + + if (ierr) { + if (verbose || veryverbose) { CCTK_VInfo (CCTK_THORNSTRING, - "Iteration %d, time %g, factor %g, throttle %g, residual %g, speed %g", - iter, (double)(currenttime-starttime), (double)factor, (double)throttle, (double)norm2, (double)speed); - } else { + "Residual calculation or boundary conditions reported error; aborting solver."); + } + Util_TableSetInt (options_table, iter, "iters"); + Util_TableDeleteKey (options_table, "error"); + goto done; + } + + // Restrict and calculate norm + ierr = 0; + BEGIN_REVERSE_REFLEVEL_LOOP(cctkGH) { + const int do_every + = ipow (mglevelfact * (maxreflevelfact/reflevelfact), + reflevelpower); + if (reflevel <= solve_level && iter % do_every == 0) { + +// cout << "CJ restrict iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl; + + if (! ierr) { + if (reflevel < solve_level) { + Restrict (cctkGH); + } + } + + BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { + if (! ierr) { + ierr = applybnds (cctkGH, options_table, userdata); + } + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; + + // Initialise norm + CCTK_REAL norm_count = 0; + CCTK_REAL norm_l2 = 0; + + BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { + if (! ierr) { + + // Calculate norm + for (int n=0; n<nvars; ++n) { + const CCTK_REAL * restrict resptr + = (const CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, res[n]); + assert (resptr); + const CCTK_REAL fac = factors[n]; + assert (cctkGH->cctk_dim == 3); + for (int k=cctkGH->cctk_nghostzones[2]; + k<cctkGH->cctk_lsh[2]-cctkGH->cctk_nghostzones[2]; + ++k) { + for (int j=cctkGH->cctk_nghostzones[1]; + j<cctkGH->cctk_lsh[1]-cctkGH->cctk_nghostzones[1]; + ++j) { + for (int i=cctkGH->cctk_nghostzones[0]; + i<cctkGH->cctk_lsh[0]-cctkGH->cctk_nghostzones[0]; + ++i) { + const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k); + ++norm_count; + // TODO: scale the norm by the resolution? + norm_l2 += pow(fac * resptr[ind], 2); + } + } + } + } // for n + + } + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; + + const int sum_handle = CCTK_ReductionArrayHandle ("sum"); + assert (sum_handle >= 0); + CCTK_REAL reduce_in[2], reduce_out[2]; + reduce_in[0] = norm_count; + reduce_in[1] = norm_l2; + ierr = CCTK_ReduceLocArrayToArray1D + (cctkGH, -1, sum_handle, reduce_in, reduce_out, 2, + CCTK_VARIABLE_REAL); + norm_counts.at(reflevel) = reduce_out[0]; + norm_l2s.at(reflevel) = reduce_out[1]; + +#warning "TODO" +#if 0 + CCTK_OutputVarAs (cctkGH, "WaveToyFO::phi", "phi-ell"); + CCTK_OutputVarAs (cctkGH, "IDSWTEsimpleFO::residual", "residual-ell"); +#endif + + } + } END_REVERSE_REFLEVEL_LOOP; + + if (ierr) { + if (verbose || veryverbose) { CCTK_VInfo (CCTK_THORNSTRING, - "Iteration %d, time %g, residual %g", - iter, (double)(currenttime-starttime), (double)norm2); + "Residual calculation or boundary conditions reported error; aborting solver."); } - if (outeveryseconds > 0) { - while (nexttime <= currenttime) nexttime += outeveryseconds; + Util_TableSetInt (options_table, iter, "iters"); + Util_TableDeleteKey (options_table, "error"); + goto done; + } + + // Save old norm + const CCTK_REAL old_norm2 = norm2; + + // Calculate new norm + CCTK_REAL norm_count = 0; + CCTK_REAL norm_l2 = 0; + for (int rl=0; rl<=solve_level; ++rl) { + norm_count += norm_counts[rl]; + norm_l2 += norm_l2s[rl]; + } + norm2 = sqrt(norm_l2 / norm_count); + + // Calculate convergence speed + speed = old_norm2 / norm2; + + // Log output + if (verbose || veryverbose) { + const time_t currenttime = time(0); + if ((iter % icoarsestep == 0 && iter >= maxiters) +#if HAVE_ISNAN + || isnan(norm2) +#endif + || norm2 <= minerror + || (iter % outevery == 0 && currenttime >= nexttime)) { + if (verbose || veryverbose) { + if (did_hit_min) { + CCTK_VInfo (CCTK_THORNSTRING, + "Convergence factor became too small: artificially kept at minfactor (may lead to instability)"); + } + if (did_hit_max) { + CCTK_VInfo (CCTK_THORNSTRING, + "Convergence factor became too large: artificially kept at maxfactor (may lead to slower convergence)"); + } + did_hit_min = false; + did_hit_max = false; + } + if (veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Iteration %d, time %g, factor %g, throttle %g, residual %g, speed %g", + iter, double(currenttime-starttime), + double(factor), double(throttle), double(norm2), + double(speed)); + } else { + CCTK_VInfo (CCTK_THORNSTRING, + "Iteration %d, time %g, residual %g", + iter, double(currenttime-starttime), double(norm2)); + } + if (outeveryseconds > 0) { + while (nexttime <= currenttime) nexttime += outeveryseconds; + } } } - } - - if (iter == maxiters) break; - + #if HAVE_ISNAN - // Exit if things go bad - if (isnan(norm2)) { - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, "Encountered NaN. Aborting solver."); + // Exit if things go bad + if (isnan(norm2)) { + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Encountered NaN. Aborting solver."); + } + break; } - break; - } #endif - - // Return when desired accuracy is reached - if (norm2 <= minerror) { - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, "Done solving."); + + // Return when desired accuracy is reached + if (norm2 <= minerror) { + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, "Done solving."); + } + Util_TableSetInt (options_table, iter, "iters"); + Util_TableSetReal (options_table, norm2, "error"); + ierr = 0; + goto done; } - Util_TableSetInt (options_table, iter, "iters"); - Util_TableSetReal (options_table, norm2, "error"); - ierr = 0; - goto done; - } + + // Exit if the maximum number of iterations has been reached + if (iter % icoarsestep == 0 && iter >= maxiters) { + break; + } + + // Adjust speed + if (speed > 1.0) { + throttle = acceleration; + } else { + throttle = deceleration; + } + factor *= throttle; + if (factor < minfactor) { + did_hit_min = true; + factor = minfactor; + } + if (factor > maxfactor) { + did_hit_max = true; + factor = maxfactor; + } + + } // for iter - // Adjust speed - if (speed > 1.0) { - throttle = acceleration; - } else { - throttle = deceleration; - } - factor *= throttle; - if (factor < minfactor) { - did_hit_min = true; - factor = minfactor; - } - if (factor > maxfactor) { - did_hit_max = true; - factor = maxfactor; + // Did not solve + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, "Did not converge."); } + Util_TableSetInt (options_table, iter, "iters"); + Util_TableSetReal (options_table, norm2, "error"); + ierr = -1; + + done: ; - // Smooth - CallLocalFunction ((cGH *)cctkGH, common::smooth); - // Apply boundary conditions - common::ierr = 0; - CallLocalFunction ((cGH *)cctkGH, common::call_applybnds); - ierr = common::ierr; - if (ierr != 0) { - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, "Boundary conditions reported error; aborting solver."); - } - Util_TableSetInt (options_table, iter, "iters"); - Util_TableDeleteKey (options_table, "error"); - goto done; - } -#if 0 - // Restrict - BEGIN_REVERSE_REFLEVEL_LOOP(cctkGH) { - BEGIN_MGLEVEL_LOOP(cctkGH) { - Restrict (cctkGH); - // TODO: do something here -// CCTK_ScheduleTraverse ("CCTK_POSTRESTRICT", (cGH *)cctkGH, CallFunction); - CallLocalFunction ((cGH *)cctkGH, common::call_applybnds); - assert (! common::ierr); - } END_MGLEVEL_LOOP; - } END_REVERSE_REFLEVEL_LOOP; -#endif + // Reset the initial time +#warning "TODO: reset the initial time a bit more intelligently" + global_time = 0; + delta_time = 1; + * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; + * const_cast<CCTK_REAL *> (& cctkGH->cctk_delta_time) = delta_time; + BEGIN_REFLEVEL_LOOP(cctkGH) { + * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; + for (int m=0; m<maps; ++m) { + vtt[m]->set_time (reflevel, mglevel, 0); + vtt[m]->set_delta (reflevel, mglevel, 1.0 / reflevelfact); + } + } END_REFLEVEL_LOOP; - } // for iter - - // Did not solve - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, "Did not converge."); - } - Util_TableSetInt (options_table, iter, "iters"); - Util_TableSetReal (options_table, norm2, "error"); - ierr = -1; - - done: - - // Restore state - if (reflevel!=saved_reflevel) { - set_reflevel ((cGH *)cctkGH, saved_reflevel); - set_mglevel ((cGH *)cctkGH, saved_mglevel); - } + } END_GLOBAL_MODE; return ierr; } |