[cig-commits] commit: Output the residual

Mercurial hg at geodynamics.org
Fri Feb 10 16:00:12 PST 2012


changeset:   13:a48ee28dbd74
user:        Walter Landry <wlandry at caltech.edu>
date:        Mon Nov 28 14:32:56 2011 -0800
files:       main.cxx
description:
Output the residual


diff -r e61eb71c2e37 -r a48ee28dbd74 main.cxx
--- a/main.cxx	Mon Nov 28 14:32:45 2011 -0800
+++ b/main.cxx	Mon Nov 28 14:32:56 2011 -0800
@@ -3,6 +3,7 @@
 
 #include "Model.hxx"
 #include "constants.hxx"
+#include "write_vtk.hxx"
 
 extern void initial(const Model &model, double zx[N+1][N], double zy[N][N+1],
                     double log_etax[N+1][N], double log_etay[N][N+1],
@@ -32,7 +33,8 @@ int main()
 
   initial(model,zx,zy,log_etax,log_etay,p,fx,fy);
 
-  return 0;
+  double Resid_p[N][N], Resid_x[N+1][N], Resid_y[N][N+1];
+
   for(int iteration=0;iteration<max_iterations;++iteration)
     {
       /* Compute the boundary jumps */
@@ -89,7 +91,6 @@ int main()
         {
 
           /* zx */
-          std::ofstream outfilex("zx");
           for(int rb=0;rb<2;++rb)
             {
               for(int j=0;j<N;++j)
@@ -206,21 +207,17 @@ int main()
 
                           double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
 
+                          Resid_x[i][j]=Rx;
                           // zx[i][j]-=theta_mom*Rx/dRx_dzx;
-
-                          outfilex << x << " "
-                                   << y << " "
-                                   << Rx << " "
-                                   << "\n";
                         }
                     }
                 }
             }
 
+          write_vtk("zx_resid",N+1,N,Resid_x);
 
           /* zy */
 
-          std::ofstream outfiley("zy");
           for(int rb=0;rb<2;++rb)
             {
               for(int j=0;j<N+1;++j)
@@ -330,20 +327,17 @@ int main()
 
                           double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
 
+                          Resid_y[i][j]=Ry;
                           // zy[i][j]-=theta_mom*Ry/dRy_dzy;
-
-                          outfiley << x << " "
-                                   << y << " "
-                                   << Ry << " "
-                                   << "\n";
                         }
                     }
                 }
             }
 
+          write_vtk("zy_resid",N,N,Resid_y);
+
           /* Pressure update */
 
-          std::ofstream outfilep("p");
           for(int j=0;j<N;++j)
             for(int i=0;i<N;++i)
               {
@@ -397,14 +391,13 @@ int main()
                   + dRc_dzyp*dRmp_dp/dRm_dz
                   + dRc_dzym*dRmm_dp/dRm_dz;
 
+                Resid_p[i][j]=Rc;
+
                 dp[i][j]=-theta_c*Rc/dRc_dp;
                 // p[i][j]+=dp[i][j];
+              }
 
-                outfilep << x << " "
-                         << y << " "
-                         << Rc << " "
-                         << "\n";
-              }
+          write_vtk("p_resid",N,N,Resid_p);
 
           /* Velocity Fix */
 
@@ -479,28 +472,4 @@ int main()
           max_p_previous=max_p;
         }
     }
-  // {
-  //   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";
-  // }
 }



More information about the CIG-COMMITS mailing list