Corrected the new working example KGMinkowski.
authorAndras Laszlo <laszlo.andras@wigner.mta.hu>
Sat, 29 Aug 2015 12:39:33 +0000 (14:39 +0200)
committerAndras Laszlo <laszlo.andras@wigner.mta.hu>
Sat, 29 Aug 2015 12:39:33 +0000 (14:39 +0200)
inputs/example/kgminkowski/README [new file with mode: 0644]
inputs/example/kgminkowski/analysis/README [new file with mode: 0644]
inputs/example/kgminkowski/analysis/src/intq.cc
inputs/example/kgminkowski/kgminkowski.in
man/gridripperapp.man
src/phys/example/kgminkowski/KGMinkowski.cxx

diff --git a/inputs/example/kgminkowski/README b/inputs/example/kgminkowski/README
new file mode 100644 (file)
index 0000000..ef66b91
--- /dev/null
@@ -0,0 +1,9 @@
+
+This example produces gridripper output using KGMinkowski PDE.
+
+In order to run it, say:
+gridripper kgminkowski.in
+
+The produced .bdata files can be analyzed using the program in the 
+analysis directory. See the README file there.
+
diff --git a/inputs/example/kgminkowski/analysis/README b/inputs/example/kgminkowski/analysis/README
new file mode 100644 (file)
index 0000000..1ed8093
--- /dev/null
@@ -0,0 +1,27 @@
+
+This analysis program reads .bdata files of KGMinkowski PDE and checks 
+energy and angular momentum conservation.
+
+In order to run, first produce .bdata file with gridripper using input 
+file in the previous directory.
+
+cd ..
+gridripper kgminkowski.in
+
+Then, come back to this directory, and run the analysis program:
+
+cd analysis
+./bin/go.sh
+
+(The ./bin/go.sh script is basically equivalent to
+  make
+  ./bin/intq inp/intq.conf
+commands.)
+
+Then, view the outputs under:
+
+out/data/intq.dat
+out/data/density.dat
+
+Recommended to view e.g. with gnuplot.
+
index a06b16c..f9db3d5 100644 (file)
@@ -39,7 +39,7 @@ int main (int argc, char *argv[])
 
     // Get the name of the configfile from the argument:
     if ( argc!=2 ) { cerr<<"Usage:\tintq <configfile>\n"; exit(1); }
-    string configfile=argv[1];
+    const string configfile=argv[1];
 
     // Set gstream debug level:
     gstream_debug(1);
@@ -51,7 +51,7 @@ int main (int argc, char *argv[])
     conf.append(configfile);
 
     // Get the input file name:
-    string inputFile=getconfig(conf, "inputFile");
+    const string inputFile=getconfig(conf, "inputFile");
 
     // Append the input file parameters to configuration:
     confin.append(inputFile);
@@ -63,8 +63,8 @@ int main (int argc, char *argv[])
        if ( inputFile.substr(0, 6)=="/grid/" ) isongrid=true;
     }
     vector<string> bdataFile;
-    string::size_type lastPerIndex=inputFile.rfind("/");
-    string bdataDir=(lastPerIndex!=string::npos ? inputFile.substr(0, lastPerIndex+1) : ".");
+    const string::size_type lastPerIndex=inputFile.rfind("/");
+    const string bdataDir=(lastPerIndex!=string::npos ? inputFile.substr(0, lastPerIndex+1) : ".");
     string bdataMask=(lastPerIndex!=string::npos ? inputFile.substr(lastPerIndex+1, inputFile.length()-(lastPerIndex+1)) : inputFile);
     if ( bdataMask.length()>=3 )
     {
@@ -78,25 +78,24 @@ int main (int argc, char *argv[])
     bdataList.clear();
 
     // Get output file name for integrated quantities:
-    string intqFile=getconfig(conf, "intqFile");
+    const string intqFile=getconfig(conf, "intqFile");
     // Get output file name for density quantities:
-    string densityFile=getconfig(conf, "densityFile");
+    const string densityFile=getconfig(conf, "densityFile");
 
     // Get log file name:
-    string logFile=getconfig(conf, "logFile");
+    const string logFile=getconfig(conf, "logFile");
     summaryinfo::logfilename(logFile);
 
     // Get EPSILON:
-    GReal_t EPSILON=getconfig(conf, "EPSILON");
+    const GReal_t EPSILON=getconfig(conf, "EPSILON");
 
     // Simaulation parameters:
-    GReal_t minr=0.0;
-    GReal_t minrIntq=minr;
-    GReal_t maxr=getconfig(confin, "parameter.maxr");
-    GReal_t discr=getconfig(confin, "maxi");
-    GReal_t maxrIntq=minrIntq+(maxr-minr)/discr*int(0.8*discr);
-    int maxOrder=getconfig(confin, "parameter.maxOrder");
-    bool isTimeReversed=getconfig(confin, "parameter.isTimeReversed");
+    const GReal_t minr=0.0;
+    const GReal_t minrIntq=minr;
+    const GReal_t maxr=getconfig(confin, "parameter.maxr");
+    const GReal_t discretization=getconfig(confin, "maxi");
+    const GReal_t maxrIntq=minrIntq+(maxr-minr)/discretization*int(0.8*discretization);
+    const int maxOrder=getconfig(confin, "parameter.maxOrder");
 
     // Loop over time series:
     ScalarFieldMP<GComplex_t> physf(maxOrder);
@@ -113,53 +112,56 @@ int main (int argc, char *argv[])
           << "Eflow("<<maxrIntq<<")\t" 
           << "Ltot\t" << "L("<<minrIntq<<","<<maxrIntq<<")\t" 
           << "Lflow("<<minrIntq<<")\t" 
-          << "Lflow("<<maxrIntq<<")\t
-          <<endl;
+          << "Lflow("<<maxrIntq<<")" 
+          << endl;
     intq.flush();
     ofstream density;
     density.open(densityFile.c_str());
     density << "# " << "t\t" <<"r\t"<< "Refin\t" << "Edens\t" << "Ldens" <<endl;
+    // Loop over bdata files:
     for ( int i=0 ; i<(int)bdataFile.size() ; ++i )
     {
        BDataReader bdataReader(bdataFile[i]);
+       // Loop over time slices:
        while ( bdataReader.read() )
        {
-           GReal_t t=bdataReader.getTime();
-           int nc=bdataReader.getNumComponents();
+           const GReal_t t=bdataReader.getTime();
+           const int nc=bdataReader.getNumComponents();
            if ( nc%2!=0 || (nc/2)%3!=0 || (maxOrder+1)*(maxOrder+1)*2*3!=nc ) { cerr<<"Inconsistent number of components!\n"; exit(1); }
-           GReal_t xMin=bdataReader.getXMin();
+           const GReal_t xMin=bdataReader.getXMin();
            const tvalarray<int> &kXArray=bdataReader.getKxArray();
-           int maxi=(int)kXArray.size()-1;
-           GReal_t dx=bdataReader.getDeltaX();
-           int xMultiplier=bdataReader.getXMultiplier();
+           const int maxi=(int)kXArray.size()-1;
+           const GReal_t dx=bdataReader.getDeltaX();
+           const int xMultiplier=bdataReader.getXMultiplier();
            const GReal_t *theData=(const GReal_t*)bdataReader.getDoubleData();
            tvector<GReal_t> totalR, partialR;
            tvector<GReal_t> totalE, partialE;
            tvector<GReal_t> totalL, partialL;
+           // Loop over space points:
            for ( int j=0 ; j<=maxi ; ++j )
            {
-               GReal_t r=xMin+kXArray[j]*dx;
-               GReal_t refinementFactor=(GReal_t)xMultiplier/(GReal_t)(j<maxi ? kXArray[j+1]-kXArray[j] : kXArray[j]-kXArray[j-1]);
+               const GReal_t r=xMin+kXArray[j]*dx;
+               const GReal_t refinementFactor=(GReal_t)xMultiplier/(GReal_t)(j<maxi ? kXArray[j+1]-kXArray[j] : kXArray[j]-kXArray[j-1]);
                GReal_t *fieldAtX=(GReal_t*)(&theData[j*nc]);
-               int nc2=(maxOrder+1)*(maxOrder+1)*2;
-               ScalarFieldMP<GComplex_t> f(&fieldAtX[nc2*0], maxOrder);
-               ScalarFieldMP<GComplex_t> ft(&fieldAtX[nc2*1], maxOrder);
-               ScalarFieldMP<GComplex_t> fr(&fieldAtX[nc2*2], maxOrder);
+               const int nc2=(maxOrder+1)*(maxOrder+1)*sizeof(GComplex_t)/sizeof(GReal_t);
+               const ScalarFieldMP<GComplex_t> f(&fieldAtX[nc2*phys::example::kgminkowski::KGMinkowski::ind_f], maxOrder);
+               const ScalarFieldMP<GComplex_t> ft(&fieldAtX[nc2*phys::example::kgminkowski::KGMinkowski::ind_ft], maxOrder);
+               const ScalarFieldMP<GComplex_t> fr(&fieldAtX[nc2*phys::example::kgminkowski::KGMinkowski::ind_fr], maxOrder);
                physf=f/r;
                physft=ft/r;
                physfr=fr/r-f/(r*r);
                totalR.push_back(r);
                if ( minrIntq-EPSILON<r && r<maxrIntq+EPSILON ) partialR.push_back(r);
                // Calculate energy density:
-               GReal_t Edensity=(r<minr+EPSILON ? 0.0 :
-                                       real(0.5*(r*r)*L2S2(physft, physft)+0.5*(r*r)*L2S2(physfr, physfr)-0.5*L2S2(physf, LaplaceS2(physf)))
-                                );
+               const GReal_t Edensity=(r<minr+EPSILON ? 0.0 :
+                                        real(0.5*(r*r)*L2S2(physft, physft)+0.5*(r*r)*L2S2(physfr, physfr)-0.5*L2S2(physf, LaplaceS2(physf)))
+                                      );
                totalE.push_back(Edensity);
                if ( minrIntq-EPSILON<r && r<maxrIntq+EPSILON ) partialE.push_back(Edensity);
                // Calculate angular momentum density:
-               GReal_t Ldensity=(r<minr+EPSILON ? 0.0 :
-                                       real((r*r)*L2S2(physft, partial_phi(physf)))
-                                );
+               const GReal_t Ldensity=(r<minr+EPSILON ? 0.0 :
+                                        real((r*r)*L2S2(physft, partial_phi(physf)))
+                                      );
                totalL.push_back(Ldensity);
                if ( minrIntq-EPSILON<r && r<maxrIntq+EPSILON ) partialL.push_back(Ldensity);
                density << t << "\t" << r << "\t" << refinementFactor << "\t" << Edensity << "\t" << Ldensity <<endl;
@@ -168,16 +170,13 @@ int main (int argc, char *argv[])
                    if ( tSequence.size()==0 ) tSequence.push_back(t);
                    else if ( ::fabs(tSequence[tSequence.size()-1]-t)>EPSILON ) tSequence.push_back(t);
                    // Calculate energy flow:
-//                 GReal_t Eflow=(r<minr+EPSILON ? 0.0 :
-//                                     -Delta*real(L2S2(physfr, physft)
-//                                     -GI*q*e*r*L2S2(physf, divideByNeumann(scaledSigma, physfr, normBoundNeumann(scaledSigma(0, 0), H2S2SigmaSigma/norm(scale)), neumannTolerance))/scale)
-//                               );
-                   GReal_t Eflow=0.0;
+                   const GReal_t Eflow=(r<minr+EPSILON ? 0.0 :
+                                         -(r*r)*real(L2S2(physfr, physft))
+                                       );
                    // Calculate angular momentum flow:
-//                 GReal_t Lflow=(r<minr+EPSILON ? 0.0 :
-//                                     -real(Delta*L2S2(physfr, partial_phi(physf))+GI*q*a*e*r*L2S2(physf, divideByNeumann(scaledSigma, mulSinThetaSqr(physfr), normBoundNeumann(scaledSigma(0, 0), H2S2SigmaSigma/norm(scale)), neumannTolerance)/scale))
-//                               );
-                   GReal_t Lflow=0.0;
+                   const GReal_t Lflow=(r<minr+EPSILON ? 0.0 :
+                                         -(r*r)*real(L2S2(physfr, partial_phi(physf)))
+                                       );
                    if ( ::fabs(r-minrIntq)<EPSILON )
                    {
                        EfluxInner.push_back(Eflow);
@@ -198,7 +197,7 @@ int main (int argc, char *argv[])
                   << o5integrator(tSequence, EfluxOuter) << "\t"
                   << o5integrator(totalR, totalL) << "\t" << o5integrator(partialR, partialL) << "\t" 
                   << o5integrator(tSequence, LfluxInner) << "\t"
-                  << o5integrator(tSequence, LfluxOuter) << "\t"
+                  << o5integrator(tSequence, LfluxOuter)
                   << endl;
            intq.flush();
        }
index a433053..e457fe4 100644 (file)
@@ -7,18 +7,18 @@ initCond=example.kgminkowski.Hunch
 parameters=t0,amplitudeF,centerF,widthF,rampF,lF:%d,mF:%d,amplitudeFt,centerFt,widthFt,rampFt,lFt:%d,mFt:%d,omega,maxOrder:%d,maxr,isTimeReversed:%d
 parameter.t0=0.0
 parameter.amplitudeF=1.0
-parameter.centerF=31.9832359852999
-parameter.widthF=34.3219506001479
-parameter.rampF=34.3219506001479
+parameter.centerF=32.0
+parameter.widthF=16.0
+parameter.rampF=16.0
 parameter.lF=2
 parameter.mF=2
 parameter.amplitudeFt=0.0
-parameter.centerFt=31.9832359852999
-parameter.widthFt=34.3219506001479
-parameter.rampFt=34.3219506001479
+parameter.centerFt=32.0
+parameter.widthFt=16.0
+parameter.rampFt=16.0
 parameter.lFt=2
 parameter.mFt=2
-parameter.omega=0.433804363739061
+parameter.omega=0.5
 parameter.maxOrder=6
 parameter.maxr=64.0
 parameter.isTimeReversed=0
index f2654dd..569880b 100644 (file)
@@ -85,11 +85,13 @@ and
 
 access path. This example solves the wave equation in 3+1 dimensional
 Minkowski spacetime for fields expanded with respect to spherical
-harmonics.
+harmonics. Also an analysis code is shipped with this example, for
+checking energy conservation. The analysis code is located under the
+above mentioned directory of input file.
 
-The binary-to-ASCII converter program
+A generic binary-to-ASCII converter program
 .B gettslice
-which is shipped with gridridripper+gridripperapp, must be started
+is shipped with gridridripper+gridripperapp, and can be started
 with proper command line arguments. For a help on these, just start
 it on the command line without arguments.
 
index 05ef983..c37255b 100644 (file)
@@ -81,7 +81,7 @@ KGMinkowski::KGMinkowski(const Parameters* p, const int maxorder, const int arra
     PDE(NUM_COMPS*arraysize, NUM_INDEP_COMPS*arraysize, 
        componentNames(maxorder), constraintNames(maxOrder), COORD_LABELS),
     maxOrder(maxorder), arraySize(arraysize),
-    I_f(ind_f*arraySize), I_ft(ind_ft*arraySize), I_fr(ind_f*arraySize),
+    I_f(ind_f*arraySize), I_ft(ind_ft*arraySize), I_fr(ind_fr*arraySize),
     df(NUM_COMPS*arraySize), fieldat(NUM_AUX_POINTS)
 {
     const string location="gridripper::phys::example::kgminkowski::KGMinkowski::KGMinkowski(const Parameters*,.....)";