[cig-commits] commit: Solve for v on the interface iteratively. Works for 10^6 viscosity variation
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:28 PST 2012
changeset: 29:7904eebd094d
user: Walter Landry <wlandry at caltech.edu>
date: Tue Feb 07 05:23:33 2012 -0800
files: compute_v_on_interface.cxx constants.hxx
description:
Solve for v on the interface iteratively. Works for 10^6 viscosity variation
diff -r 5366956fdeb5 -r 7904eebd094d compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Tue Feb 07 03:31:22 2012 -0800
+++ b/compute_v_on_interface.cxx Tue Feb 07 05:23:33 2012 -0800
@@ -1,3 +1,7 @@
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+
#include <gsl/gsl_spline.h>
#include "constants.hxx"
@@ -26,6 +30,10 @@ void compute_v_on_interface(double zx[N+
for(int i=0;i<N+1;++i)
y_y[i]=h*i;
+ const double a1(h), a2(2*a1), b1(a1), b2(a2);
+ double ux1_p[N], ux2_p[N], ux1_m[N], ux2_m[N],
+ uy1_p[N+1], uy2_p[N+1], uy1_m[N+1], uy2_m[N+1];
+
{
int i(middle/h);
double dx=middle/h-i;
@@ -43,6 +51,16 @@ void compute_v_on_interface(double zx[N+
double high=interpolate_v(zx,log_etax,i+2,j,h,(-2+dx)*h);
vx[j]=low*(1-dx) + high*dx;
+
+ ux1_p[j]=zx[i+1][j]*exp(-log_etax[i+1][j])*(1-dx)
+ + zx[i+2][j]*exp(-log_etax[i+2][j])*dx;
+ ux2_p[j]=zx[i+2][j]*exp(-log_etax[i+2][j])*(1-dx)
+ + zx[i+3][j]*exp(-log_etax[i+3][j])*dx;
+
+ ux1_m[j]=zx[i-1][j]*exp(-log_etax[i-1][j])*(1-dx)
+ + zx[i][j]*exp(-log_etax[i][j])*dx;
+ ux2_m[j]=zx[i-2][j]*exp(-log_etax[i-2][j])*(1-dx)
+ + zx[i-1][j]*exp(-log_etax[i-1][j])*dx;
}
}
{
@@ -62,38 +80,116 @@ void compute_v_on_interface(double zx[N+
double high=interpolate_v(zy,log_etay,i+2,j,h,h*(-2+dx));
vy[j]=low*(1-dx) + high*dx;
+
+ uy1_p[j]=zy[i+1][j]*exp(-log_etay[i+1][j])*(1-dx)
+ + zy[i+2][j]*exp(-log_etay[i+2][j])*dx;
+ uy2_p[j]=zy[i+2][j]*exp(-log_etay[i+2][j])*(1-dx)
+ + zy[i+3][j]*exp(-log_etay[i+3][j])*dx;
+
+ uy1_m[j]=zy[i-1][j]*exp(-log_etay[i-1][j])*(1-dx)
+ + zy[i][j]*exp(-log_etay[i][j])*dx;
+ uy2_m[j]=zy[i-2][j]*exp(-log_etay[i-2][j])*(1-dx)
+ + zy[i-1][j]*exp(-log_etay[i-1][j])*dx;
}
}
- /* Compute dvx_dy, dvx_dyy */
- gsl_interp_accel *vx_accel=gsl_interp_accel_alloc ();
- gsl_spline *vx_spline=gsl_spline_alloc(gsl_interp_cspline, N);
- gsl_spline_init(vx_spline, y_x, vx, N);
+ /* Iterate to find v on the boundary */
- dvx_y[0]=0;
- dvx_y[N]=0;
- for(int i=1;i<N;++i)
- dvx_y[i]=gsl_spline_eval_deriv(vx_spline,y_y[i],vx_accel);
+ const double relax=0.25;
+ const double tolerance=1e-6;
+ double error=0;
+ // for(int ii=0;ii<100;++ii)
+ do
+ {
+ error=0;
+ /* Compute the derivatives from the current guess */
+ /* dvx_dy, dvx_dyy */
+ gsl_interp_accel *vx_accel=gsl_interp_accel_alloc ();
+ gsl_spline *vx_spline=gsl_spline_alloc(gsl_interp_cspline, N);
+ gsl_spline_init(vx_spline, y_x, vx, N);
- for(int i=0;i<N;++i)
- dvx_yy[i]=gsl_spline_eval_deriv2(vx_spline,y_x[i],vx_accel);
+ dvx_y[0]=0;
+ dvx_y[N]=0;
+ for(int i=1;i<N;++i)
+ dvx_y[i]=gsl_spline_eval_deriv(vx_spline,y_y[i],vx_accel);
- gsl_spline_free(vx_spline);
- gsl_interp_accel_free(vx_accel);
+ for(int i=0;i<N;++i)
+ dvx_yy[i]=gsl_spline_eval_deriv2(vx_spline,y_x[i],vx_accel);
- /* Compute dvy_dy, dvy_dyy */
- gsl_interp_accel *vy_accel=gsl_interp_accel_alloc ();
- gsl_spline *vy_spline=gsl_spline_alloc(gsl_interp_cspline, N+1);
- gsl_spline_init(vy_spline, y_y, vy, N+1);
+ gsl_spline_free(vx_spline);
+ gsl_interp_accel_free(vx_accel);
- dvy_y[0]=0;
- dvy_y[N]=0;
- for(int i=1;i<N+1;++i)
- dvy_y[i]=gsl_spline_eval_deriv(vy_spline,y_x[i],vy_accel);
+ /* dvy_dy, dvy_dyy */
+ gsl_interp_accel *vy_accel=gsl_interp_accel_alloc ();
+ gsl_spline *vy_spline=gsl_spline_alloc(gsl_interp_cspline, N+1);
+ gsl_spline_init(vy_spline, y_y, vy, N+1);
- for(int i=0;i<N;++i)
- dvy_yy[i]=gsl_spline_eval_deriv2(vy_spline,y_y[i],vy_accel);
+ dvy_y[0]=0;
+ dvy_y[N]=0;
+ for(int i=1;i<N+1;++i)
+ dvy_y[i]=gsl_spline_eval_deriv(vy_spline,y_x[i],vy_accel);
- gsl_spline_free(vy_spline);
- gsl_interp_accel_free(vy_accel);
+ for(int i=0;i<N;++i)
+ dvy_yy[i]=gsl_spline_eval_deriv2(vy_spline,y_y[i],vy_accel);
+
+ gsl_spline_free(vy_spline);
+ gsl_interp_accel_free(vy_accel);
+
+ /* Update the current guess based on the derivatives */
+ for(int i=0;i<N;++i)
+ {
+ double zx_jump=-eta_jump*dvy_y[i];
+ double temp=
+ (min_eta*(a1*a2/(b2-b1))*(b2*b2*ux1_m[i] - b1*b1*ux2_m[i])
+ + max_eta*(b1*b2/(a2-a1))*(a2*a2*ux1_p[i] - a1*a1*ux2_p[i]))
+ /(min_eta*a1*a2*(b1+b2) + max_eta*b1*b2*(a1+a2))
+ - a1*a2*b1*b2*zx_jump/(min_eta*a1*a2*(b1+b2)
+ + max_eta*b1*b2*(a1+a2));
+ if(i!=0 && i!=N-1)
+ {
+ error=
+ std::max(error,std::fabs((vx[i]-temp)
+ /(std::fabs(vx[i])+std::fabs(temp))));
+ if(i==100)
+ std::cout << "vx error "
+ << i << " "
+ << vx[i] << " "
+ << temp << " "
+ << error << " "
+ << "\n";
+ }
+ vx[i]=(1-relax)*vx[i] + relax*temp;
+ }
+
+ for(int i=0;i<N+1;++i)
+ {
+ double zy_jump=-eta_jump*dvx_y[i];
+ double temp=
+ (min_eta*(a1*a2/(b2-b1))*(b2*b2*uy1_m[i] - b1*b1*uy2_m[i])
+ + max_eta*(b1*b2/(a2-a1))*(a2*a2*uy1_p[i] - a1*a1*uy2_p[i]))
+ /(min_eta*a1*a2*(b1+b2) + max_eta*b1*b2*(a1+a2))
+ - a1*a2*b1*b2*zy_jump/(min_eta*a1*a2*(b1+b2)
+ + max_eta*b1*b2*(a1+a2));
+ if(i!=0 && i!=N)
+ {
+ error=
+ std::max(error,std::fabs((vy[i]-temp)
+ /(std::fabs(vy[i])+std::fabs(temp))));
+ if(i==100)
+ std::cout << "vy error "
+ << i << " "
+ << vy[i] << " "
+ << temp << " "
+ << std::fabs((vy[i]-temp)
+ /(std::fabs(vy[i])+std::fabs(temp))) << " "
+ << error << " "
+ << "\n";
+ }
+ vy[i]=(1-relax)*vy[i] + relax*temp;
+ }
+ std::cout << "error "
+ << error << " "
+ << "\n";
+ }
+ while (error>tolerance);
}
diff -r 5366956fdeb5 -r 7904eebd094d constants.hxx
--- a/constants.hxx Tue Feb 07 03:31:22 2012 -0800
+++ b/constants.hxx Tue Feb 07 05:23:33 2012 -0800
@@ -1,6 +1,6 @@ const int N(1024);
-const int N(1024);
+const int N(128);
const double min_eta=1;
-const double max_eta=1e4;
+const double max_eta=1e6;
const double eta_jump=max_eta-min_eta;
#include <cmath>
const double log_max_eta=std::log(max_eta);
More information about the CIG-COMMITS
mailing list