[cig-commits] commit: FD Derivatives for zx seem to be correct
Mercurial
hg at geodynamics.org
Fri Feb 10 15:59:57 PST 2012
changeset: 0:52dc23ba6718
user: Walter Landry <wlandry at caltech.edu>
date: Wed Aug 03 23:29:26 2011 -0700
files: main.cxx
description:
FD Derivatives for zx seem to be correct
diff -r 000000000000 -r 52dc23ba6718 main.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/main.cxx Wed Aug 03 23:29:26 2011 -0700
@@ -0,0 +1,209 @@
+#include <cmath>
+#include <iostream>
+
+enum class Model {solCx, solKz, sinker};
+
+int main()
+{
+ const int N(64);
+
+ double zx[N+1][N], zy[N][N+1], log_etax[N+1][N], log_etay[N][N+1],
+ p[N][N], fx[N+1][N], fy[N][N+1];
+
+
+ /* Initial conditions */
+
+ const double max_eta=2;
+ const double log_max_eta=std::log(max_eta);
+ const double middle=0.4;
+ const int n_sweeps=1;
+ const double h(1.0/N);
+ const double pi(atan(1.0)*4);
+
+ Model model(Model::solCx);
+
+ for(int i=0;i<N+1;++i)
+ for(int j=0;j<N+1;++j)
+ {
+ if(j<N)
+ {
+ double x(i*h), y((j+0.5)*h);
+ zx[i][j]=1.2 + 3.4*x + 4.5*y;
+ // zx[i][j]=0;
+ fx[i][j]=0;
+ if(model==Model::solCx)
+ {
+ log_etax[i][j]=y*log_max_eta;
+ }
+ else if(model==Model::solKz)
+ {
+ if(x>middle)
+ {
+ log_etax[i][j]=log_max_eta;
+ }
+ else
+ {
+ log_etax[i][j]=0;
+ }
+ }
+ else if(model==Model::sinker)
+ {
+ if(x>.4 && x < .6 && y>.2 && y<.4)
+ {
+ log_etax[i][j]=log_max_eta;
+ }
+ else
+ {
+ log_etax[i][j]=0;
+ }
+ }
+ }
+
+ if(i<N)
+ {
+ double x((i+0.5)*h), y(j*h);
+ 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)
+ {
+ fy[i][j]=sin(pi*y)*cos(pi*x);
+ if(model==Model::solCx)
+ {
+ log_etay[i][j]=y*log_max_eta;
+ }
+ else
+ {
+ if(x>middle)
+ {
+ log_etay[i][j]=log_max_eta;
+ }
+ else
+ {
+ log_etay[i][j]=0;
+ }
+ }
+ }
+ else if(model==Model::sinker)
+ {
+ if(x>.4 && x < .6 && y>.2 && y<.4)
+ {
+ log_etay[i][j]=log_max_eta;
+ fy[i][j]=1;
+ }
+ else
+ {
+ log_etay[i][j]=0;
+ fy[i][j]=0;
+ }
+ }
+ }
+
+ if(i<N && j<N)
+ p[i][j]=0;
+ }
+
+ for(int sweep=0;sweep<n_sweeps;++sweep)
+ {
+ for(int rb=0;rb<2;++rb)
+ {
+ for(int j=0;j<N;++j)
+ {
+ int i_min=(j+rb)%2;
+ for(int i=i_min;i<N+1;i+=2)
+ {
+ if(i==0 || i==N)
+ {
+ zx[i][j]=0;
+ }
+ else
+ {
+ /* Do the finite difference parts */
+
+ /* Derivatives of z */
+ double dzx_xx=(zx[i-1][j] - 2*zx[i][j] + zx[i+1][j])/(h*h);
+
+ double dzy_xy=((zy[i][j+1] - zy[i-1][j+1])
+ - (zy[i][j] - zy[i-1][j]))/(h*h);
+
+ double dzx_x=(zx[i+1][j] - zx[i-1][j])/(2*h);
+ double dzy_y=(zy[i][j+1] - zy[i][j]
+ + zy[i-1][j+1] - zy[i-1][j])/(2*h);
+
+ double dzx_y, dzx_yy;
+ if(j==0)
+ {
+ dzx_yy=(-zx[i][j] + zx[i][j+1])/(h*h);
+ dzx_y=(zx[i][j+1] - zx[i][j])/(2*h);
+ }
+ else if(j==N-1)
+ {
+ dzx_yy=(zx[i][j-1] - zx[i][j])/(h*h);
+ dzx_y=(zx[i][j] - zx[i][j-1])/(2*h);
+ }
+ else
+ {
+ dzx_yy=(zx[i][j-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
+ 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. */
+
+ double dp_x=(p[i][j]-p[i-1][j])/h;
+
+ double dlog_etaxx=
+ ((log_etay[i][j+1] + log_etay[i][j])/2
+ - 2*log_etax[i][j]
+ + (log_etay[i-1][j+1] + log_etay[i-1][j])/2)/(h*h/4);
+ double dlog_etax=
+ ((log_etay[i][j+1] + log_etay[i][j])/2
+ - (log_etay[i-1][j+1] + log_etay[i-1][j])/2)/h;
+
+ double dlog_etayy=
+ ((log_etay[i][j+1] + log_etay[i-1][j+1])/2
+ - 2*log_etax[i][j]
+ + (log_etay[i][j] + log_etay[i-1][j])/2)/(h*h/4);
+
+ double dlog_etay=
+ ((log_etay[i][j+1] + log_etay[i-1][j+1])
+ - (log_etay[i][j] + log_etay[i-1][j]))/(2*h);
+
+ double dlog_etaxy=
+ ((log_etay[i][j+1] - log_etay[i-1][j+1])
+ - (log_etay[i][j] - log_etay[i-1][j]))/(h*h);
+ /* Now do the jump conditions */
+
+ // const double x(i*h), y((j+0.5)*h);
+ // std::cout << i << " "
+ // << j << " "
+ // << x << " "
+ // << y << " "
+ // // << h << " "
+ // << zx[i][j] << " "
+ // << dzx_x << " "
+ // << dzx_xx << " "
+ // << dzx_y << " "
+ // << dzx_yy << " "
+ // << dzy_y << " "
+ // << dzy_xy << " "
+
+ // // << log_etax[i][j] << " "
+ // // << log_max_eta << " "
+ // // << dlog_etax << " "
+ // // << dlog_etay << " "
+ // // << dlog_etaxx << " "
+ // // << dlog_etayy << " "
+ // // << dlog_etaxy << " "
+ // << "\n";
+
+ }
+ }
+ }
+ }
+ }
+}
More information about the CIG-COMMITS
mailing list