aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorschnetter <>2004-01-25 12:57:00 +0000
committerschnetter <>2004-01-25 12:57:00 +0000
commit55f76be95c2272b35b6941d6ec38a77bfb23a101 (patch)
treef1c450c6d49389fb5ef2c94341138c003bd594ef
parent035e2124aec729716b9c7db5ef0e1a92f1ea0d6a (diff)
Import the recently announced changes:
Import the recently announced changes: 1. Carpet has now an infrastructure for multiple maps (aka "grid patches"). Instead of a single grid hierarchy there can now be several. This is largely untested, because the remainder of Cactus cannot handle multiple coordinate systems. 2. The order in which the schedule bins are called has changed. As Ian Hawke pointed out, the previous order during time evolution was inconsistent. The initial data ordering did not allow for recovering and was not usable for progressively solving elliptic equations for initial data. 3. Carpet now supports convergence levels. The convergence level specifies by how many factors of two the resolution in the parameter file should be coarsened (or refined, if negative). This should make convergence tests and test runs much easier. It is, in principle, also possible to run several convergence levels at once. This has not been tested because the remainder of Cactus cannot handle multiple resolutions. This will be necessary for a multigrid solver, and also for having a shadow hierarchy to determine where to refine adaptively. 4. Carpet works together with the new CoordBase domain specification parameters. Without these, using convergence levels will lead to very strange results. 5. The "modes" have changed. There are now: meta mode: the whole simulation global mode: one convergence level level mode: one refinement level singlemap mode: one map on one refinement level local mode: as previously The whole mode handling has been cleaned up. 6. The regridding thorn has been cleaned up. 7. The kind of prolongation stencil is now determined in Carpet, i.e. at a fairly hight level, instead of in CarpetLib. 8. The low-order prolongation operators have been made much more efficient (as have previously the higher-order ones). 9. Assorted smaller changes. For Carpet users, there should be no major incompatibilities. The major improvements are 3 and 4 combined. Here is an example: CoordBase::domainsize = extent CoordBase::spacing = gridspacing CoordBase::zero_origin_x = yes CoordBase::zero_origin_y = yes CoordBase::zero_origin_z = yes CoordBase::xextent = 20.0 CoordBase::yextent = 20.0 CoordBase::zextent = 20.0 CoordBase::dx = 1.0 CoordBase::dy = 1.0 CoordBase::dz = 1.0 CoordBase::boundary_shiftout_x_lower = 1 CoordBase::boundary_shiftout_y_lower = 1 CoordBase::boundary_shiftout_z_lower = 1 Carpet::domain_from_coordbase = yes Carpet::convergence_level = 0 grid::type = coordbase grid::domain = octant grid::avoid_origin = no This gives you a grid that extends from the origin ("zero_origin") up to 20.0 with a grid spacing of 1.0. Symmetry zones and boundary zones are added automatically. The "shiftout" says that there is no boundary point on the origin. The staggering parameters (not shown) default to "no". In order to change the resolution, only the convergence level has to be adjusted. Note that the old way of specifying the domain extent still works. For Carpet developers, one major change is the new mode handling. As described in 5, the looping macros (that loop over all refinement levels, or all components) have changed. darcs-hash:20040125125727-07bb3-3368611314b2dcb8c8ae58ab3f501b683d7edb8f.gz
Diffstat
-rw-r--r--Carpet/Carpet/README5
-rw-r--r--Carpet/Carpet/src/Shutdown.cc47
-rw-r--r--Carpet/Carpet/src/carpet_public.hh299
-rw-r--r--Carpet/Carpet/src/make.code.defn3
-rw-r--r--Carpet/CarpetIOASCII/README3
-rw-r--r--Carpet/CarpetInterp/README3
-rw-r--r--Carpet/CarpetLib/README3
-rw-r--r--Carpet/CarpetRegrid/schedule.ccl9
-rw-r--r--Carpet/CarpetSlab/README3
-rw-r--r--Carpet/CarpetSlab/src/slab.h22
-rw-r--r--Carpet/CarpetTest/src/carpettest_check_arguments.F7716
-rw-r--r--CarpetAttic/CarpetIOFlexIO/README4
-rw-r--r--CarpetDev/CarpetCG/src/CG.cc872
-rw-r--r--CarpetDev/CarpetJacobi/param.ccl14
-rw-r--r--CarpetDev/CarpetJacobi/src/Jacobi.cc698
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;
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-18 07:38:22 +0000