[cig-commits] commit: add 2D elastic smoother
Mercurial
hg at geodynamics.org
Thu Jun 7 13:35:35 PDT 2012
changeset: 256:96b2e060a8ce
user: Sylvain Barbot <sbarbot at caltech.edu>
date: Tue Apr 10 12:28:04 2012 -0700
files: input/sinker.input src/StokesFACOps.h src/StokesFACOps/residual_2D.C src/StokesFACOps/smooth_Tackley_2D.C src/StokesFACOps/smooth_V_2D.C src/dRm_dv.h
description:
add 2D elastic smoother
diff -r 567a621add59 -r 96b2e060a8ce input/sinker.input
--- a/input/sinker.input Tue Apr 10 11:15:07 2012 -0700
+++ b/input/sinker.input Tue Apr 10 12:28:04 2012 -0700
@@ -31,20 +31,12 @@ FACStokes {
adaption_threshold = 1.0e-3
min_full_refinement_level = 3
- lambda_ijk= 11, 11
+ lambda_ijk= 2, 2
lambda_coord_min= -0.001, -0.001
lambda_coord_max= 1.001, 1.001
- lambda_data= 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
+ lambda_data=
+ 1, 1,
+ 1, 1
mu_ijk= 3, 3
mu_coord_min= -0.001, -0.001
diff -r 567a621add59 -r 96b2e060a8ce src/StokesFACOps.h
--- a/src/StokesFACOps.h Tue Apr 10 11:15:07 2012 -0700
+++ b/src/StokesFACOps.h Tue Apr 10 12:28:04 2012 -0700
@@ -499,7 +499,6 @@ void smooth_V_2D(const int &axis,
const pdat::CellIndex ¢er,
const hier::Index &ip,
const hier::Index &jp,
- pdat::CellData<double> &p,
pdat::SideData<double> &v,
pdat::SideData<double> &v_rhs,
double &maxres,
@@ -524,10 +523,17 @@ void smooth_V_2D(const int &axis,
const double &dx,
const double &dy)
{
+ return
+ edge_moduli(edge+jp,1)*(v(x+jp)-v(x ))/(dy*dy)
+ - edge_moduli(edge ,1)*(v(x )-v(x-jp))/(dy*dy)
+ + (edge_moduli(edge+jp,0)+edge_moduli(edge+jp,1))*(v(y+jp)-v(y+jp-ip))/(dx*dy)
+ - (edge_moduli(edge ,0)+edge_moduli(edge ,1))*(v(y )-v(y-ip ))/(dx*dy);
+ /*
return ((v(x+jp)-v(x))/(dy*dy)
- + (v(y+jp)-v(y+jp-ip))/(dx*dy)) * edge_moduli(edge+jp)
- - ((v(x)-v(x-jp))/(dy*dy)
- + (v(y)-v(y-ip))/(dx*dy)) * edge_moduli(edge);
+ + (v(y+jp)-v(y+jp-ip))/(dx*dy)) * edge_moduli(edge+jp,1)
+ - ((v(x)-v(x-jp))/(dy*dy)
+ + (v(y)-v(y-ip))/(dx*dy)) * edge_moduli(edge,1);
+ */
}
/* The action of the velocity operator. It is written from the
@@ -535,7 +541,6 @@ void smooth_V_2D(const int &axis,
etc. to get vy. */
double v_operator_2D(pdat::SideData<double> &v,
- pdat::CellData<double> &p,
pdat::CellData<double> &cell_moduli,
pdat::NodeData<double> &edge_moduli,
const pdat::CellIndex ¢er,
@@ -547,12 +552,11 @@ void smooth_V_2D(const int &axis,
const double &dx,
const double &dy)
{
- double dtau_xx_dx=2*((v(x+ip)-v(x))*cell_moduli(center)
- - (v(x)-v(x-ip))*cell_moduli(center-ip))/(dx*dx);
+ double dtau_xx_dx=( (v(x+ip)-v(x ))*(cell_moduli(center ,0)+2*cell_moduli(center ,1))
+ -(v(x )-v(x-ip))*(cell_moduli(center-ip,0)+2*cell_moduli(center-ip,1)))/(dx*dx);
double dtau_xy_dy=dtau_mixed(v,edge_moduli,x,y,edge,ip,jp,dx,dy);
- double dp_dx=(p(center)-p(center-ip))/dx;
- return dtau_xx_dx + dtau_xy_dy - dp_dx;
+ return dtau_xx_dx + dtau_xy_dy;
}
double v_operator_3D(pdat::SideData<double> &v,
diff -r 567a621add59 -r 96b2e060a8ce src/StokesFACOps/residual_2D.C
--- a/src/StokesFACOps/residual_2D.C Tue Apr 10 11:15:07 2012 -0700
+++ b/src/StokesFACOps/residual_2D.C Tue Apr 10 12:28:04 2012 -0700
@@ -57,7 +57,7 @@ void SAMRAI::solv::StokesFACOps::residua
else
{
v_resid(x)=v_rhs(x)
- - v_operator_2D(v,p,cell_moduli,edge_moduli,center,
+ - v_operator_2D(v,cell_moduli,edge_moduli,center,
edge,x,y,ip,jp,dx,dy);
}
}
@@ -74,7 +74,7 @@ void SAMRAI::solv::StokesFACOps::residua
else
{
v_resid(y)=v_rhs(y)
- - v_operator_2D(v,p,cell_moduli,edge_moduli,center,
+ - v_operator_2D(v,cell_moduli,edge_moduli,center,
edge,y,x,jp,ip,dy,dx);
}
}
diff -r 567a621add59 -r 96b2e060a8ce src/StokesFACOps/smooth_Tackley_2D.C
--- a/src/StokesFACOps/smooth_Tackley_2D.C Tue Apr 10 11:15:07 2012 -0700
+++ b/src/StokesFACOps/smooth_Tackley_2D.C Tue Apr 10 12:28:04 2012 -0700
@@ -113,7 +113,7 @@ void SAMRAI::solv::StokesFACOps::smooth_
/* Update v */
smooth_V_2D(0,pbox,geom,center,ip,jp,
- p,v,v_rhs,maxres,dx,dy,cell_moduli,
+ v,v_rhs,maxres,dx,dy,cell_moduli,
edge_moduli,theta_momentum);
}
}
@@ -167,7 +167,7 @@ void SAMRAI::solv::StokesFACOps::smooth_
/* Update v */
smooth_V_2D(1,pbox,geom,center,jp,ip,
- p,v,v_rhs,maxres,dy,dx,cell_moduli,
+ v,v_rhs,maxres,dy,dx,cell_moduli,
edge_moduli,theta_momentum);
}
}
diff -r 567a621add59 -r 96b2e060a8ce src/StokesFACOps/smooth_V_2D.C
--- a/src/StokesFACOps/smooth_V_2D.C Tue Apr 10 11:15:07 2012 -0700
+++ b/src/StokesFACOps/smooth_V_2D.C Tue Apr 10 12:28:04 2012 -0700
@@ -14,7 +14,6 @@ void SAMRAI::solv::StokesFACOps::smooth_
const pdat::CellIndex ¢er,
const hier::Index &ip,
const hier::Index &jp,
- pdat::CellData<double> &p,
pdat::SideData<double> &v,
pdat::SideData<double> &v_rhs,
double &maxres,
@@ -57,7 +56,7 @@ void SAMRAI::solv::StokesFACOps::smooth_
edge+jp,edge,dx,dy);
double delta_Rx=v_rhs(x)
- - v_operator_2D(v,p,cell_moduli,edge_moduli,center,
+ - v_operator_2D(v,cell_moduli,edge_moduli,center,
edge,x,y,ip,jp,dx,dy);
/* No scaling here, though there should be. */
diff -r 567a621add59 -r 96b2e060a8ce src/dRm_dv.h
--- a/src/dRm_dv.h Tue Apr 10 11:15:07 2012 -0700
+++ b/src/dRm_dv.h Tue Apr 10 12:28:04 2012 -0700
@@ -18,8 +18,9 @@ double dRm_dv_2D(SAMRAI::pdat::CellData<
const double &dx,
const double &dy)
{
- return -2*(cell_moduli(center) + cell_moduli(left))/(dx*dx)
- - (edge_moduli(up_e) + edge_moduli(center_e))/(dy*dy);
+ return -( (cell_moduli(center,0) + cell_moduli(left,0))
+ +2*(cell_moduli(center,1) + cell_moduli(left,1)))/(dx*dx)
+ -(edge_moduli(up_e,1) + edge_moduli(center_e,1))/(dy*dy);
}
#endif
More information about the CIG-COMMITS
mailing list