[cig-commits] commit: Fix a mixup between solKz and solCx. Also, it does the full solKz benchmark now

[Mercurial]

changeset:   7:aa1b9fba90b4
user:        Walter Landry <wlandry at caltech.edu>
date:        Thu Aug 04 06:54:40 2011 -0700
files:       main.cxx
description:
Fix a mixup between solKz and solCx.  Also, it does the full solKz benchmark now


diff -r da6fd69b973b -r aa1b9fba90b4 main.cxx
--- a/main.cxx	Thu Aug 04 02:33:16 2011 -0700
+++ b/main.cxx	Thu Aug 04 06:54:40 2011 -0700
@@ -1,11 +1,18 @@
 #include <cmath>
 #include <iostream>
+#include <fstream>
 
-enum class Model {solCx, solKz, sinker};
+enum class Model {solKz, solCx, sinker};
 
 int main()
 {
   const int N(64);
+
+  /* Eta is collocated with z.  Maybe eta should be cell & vertex
+     centered?  It would simplify most of these differences, but it
+     would make dlog_etaxy for the momentum residual and dlog_etax for
+     the continuity residual more messy.  Also, with it collocated
+     with z, it is easy to compute v=z/eta. */
 
   double zx[N+1][N], zy[N][N+1], log_etax[N+1][N], log_etay[N][N+1],
     p[N][N], dp[N][N], fx[N+1][N], fy[N][N+1];
@@ -17,13 +24,13 @@ int main()
   const double log_max_eta=std::log(max_eta);
   const double middle=0.4;
   // const int n_sweeps=1;
-  const int n_sweeps=100000;
+  const int n_sweeps=10000;
   const double h(1.0/N);
   const double pi(atan(1.0)*4);
   const double theta_mom=1.0;
   const double theta_c=1.0;
 
-  Model model(Model::solCx);
+  Model model(Model::solKz);
 
   for(int i=0;i<N+1;++i)
     for(int j=0;j<N+1;++j)
@@ -32,19 +39,19 @@ int main()
           {
             double x(i*h), y((j+0.5)*h);
 
-            if(i==0 || i==N)
-              zx[i][j]=0;
-            else
-              zx[i][j]=1.2 + 3.4*x + 4.5*y;
+            // if(i==0 || i==N)
+            //   zx[i][j]=0;
+            // else
+            //   zx[i][j]=1.2 + 3.4*x + 4.5*y;
 
-            // zx[i][j]=0;
+            zx[i][j]=0;
 
             fx[i][j]=0;
-            if(model==Model::solCx)
+            if(model==Model::solKz)
               {
                 log_etax[i][j]=y*log_max_eta;
               }
-            else if(model==Model::solKz)
+            else if(model==Model::solCx)
               {
                 if(x>middle)
                   {
@@ -72,17 +79,17 @@ int main()
           {
             double x((i+0.5)*h), y(j*h);
 
-            if(j==0 || j==N)
-              zy[i][j]=0;
-            else
-              zy[i][j]=5.6 + 7.8*x + 9.0*y + 2.1*x*y;
-            // zy[i][j]=0;
+            // if(j==0 || j==N)
+            //   zy[i][j]=0;
+            // else
+            //   zy[i][j]=5.6 + 7.8*x + 9.0*y + 2.1*x*y;
+            zy[i][j]=0;
 
-            if(model==Model::solCx || model==Model::solKz)
+            if(model==Model::solKz || model==Model::solCx)
               {
-                fy[i][j]=0;
-                // fy[i][j]=sin(pi*y)*cos(pi*x);
-                if(model==Model::solCx)
+                // fy[i][j]=0;
+                fy[i][j]=sin(pi*y)*cos(pi*x);
+                if(model==Model::solKz)
                   {
                     log_etay[i][j]=y*log_max_eta;
                   }
@@ -122,6 +129,7 @@ int main()
 
   for(int sweep=0;sweep<n_sweeps;++sweep)
     {
+
       /* zx */
       for(int rb=0;rb<2;++rb)
         {
@@ -165,12 +173,7 @@ int main()
                           dzx_y=(zx[i][j+1] - zx[i][j-1])/(2*h);
                         }
 
-                      /* Derivatives of p and eta.  Eta is collocated
-                         with z.  Maybe eta should be cell & vertex
-                         centered?  It would simplify most of these
-                         differences, but it would make dlog_etaxy
-                         more messy.  Also, with it collocated with z,
-                         it is easy to compute v=z/eta.  */
+                      /* Derivatives of p and eta.  */
 
                       double dp_x=(p[i][j]-p[i-1][j])/h;
 
@@ -262,12 +265,7 @@ int main()
                           dzy_x=(zy[i+1][j] - zy[i-1][j])/(2*h);
                         }
 
-                      /* Derivatives of p and eta.  Eta is collocated
-                         with z.  Maxbe eta should be cell & vertey
-                         centered?  It would simplifx most of these
-                         differences, but it would make dlog_etayx
-                         more messx.  Also, with it collocated with z,
-                         it is easx to compute v=z/eta.  */
+                      /* Derivatives of p and eta. */
 
                       double dp_y=(p[i][j]-p[i][j-1])/h;
 
@@ -413,4 +411,29 @@ int main()
                   << "\n";
 
     }
+
+  {
+    std::ofstream outfile("zx");
+    for(int j=0;j<N;++j)
+      for(int i=0;i<N+1;++i)
+        outfile << i << " " << j << " "
+                << zx[i][j] << " "
+                << "\n";
+  }
+  {
+    std::ofstream outfile("zy");
+    for(int j=0;j<N+1;++j)
+      for(int i=0;i<N;++i)
+        outfile << i << " " << j << " "
+                << zy[i][j] << " "
+                << "\n";
+  }
+  {
+    std::ofstream outfile("p");
+    for(int j=0;j<N;++j)
+      for(int i=0;i<N;++i)
+        outfile << i << " " << j << " "
+                << p[i][j] << " "
+                << "\n";
+  }
 }