[cig-commits] r5980 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Feb 7 16:15:04 PST 2007
Author: tan2
Date: 2007-02-07 16:15:03 -0800 (Wed, 07 Feb 2007)
New Revision: 5980
Modified:
mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
Log:
Spherical version of the temperature solver
Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-02-08 00:14:46 UTC (rev 5979)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-02-08 00:15:03 UTC (rev 5980)
@@ -311,7 +311,7 @@
const int dims=E->mesh.nsd;
const int dofs=E->mesh.dof;
const int ends=enodes[dims];
- const int sphere_key = 0;
+ const int sphere_key = 1;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.nno;i++)
@@ -322,10 +322,14 @@
velo_from_element(E,VV,m,el,sphere_key);
- get_global_shape_fn(E,el,&GN,&GNx,&dOmega,0,sphere_key,rtf,E->mesh.levmax,m);
+ get_global_shape_fn(E, el, &GN, &GNx, &dOmega, 0,
+ sphere_key, rtf, E->mesh.levmax, m);
- pg_shape_fn(E,el,&PG,&GNx,VV,diff,m);
- element_residual(E,el,PG,GNx,dOmega,VV,T,Tdot,Q0,Eres,diff,E->sphere.cap[m].TB,FLAGS,m);
+ pg_shape_fn(E, el, &PG, &GNx, VV,
+ rtf, diff, m);
+ element_residual(E, el, PG, GNx, dOmega, VV, T, Tdot,
+ Q0, Eres, rtf, diff, E->sphere.cap[m].TB,
+ FLAGS, m);
for(a=1;a<=ends;a++) {
a1 = E->ien[m][el].node[a];
@@ -353,31 +357,25 @@
Petrov-Galerkin shape functions for a given element
=================================================== */
-void pg_shape_fn(E,el,PG,GNx,VV,diffusion,m)
+void pg_shape_fn(E,el,PG,GNx,VV,rtf,diffusion,m)
struct All_variables *E;
int el,m;
struct Shape_function *PG;
struct Shape_function_dx *GNx;
float VV[4][9];
+ double rtf[4][9];
double diffusion;
{
- int i,j,node;
+ int i,j;
int *ienm;
double uc1,uc2,uc3;
- double u1,u2,u3;
- double aa,bb,cc,uxse,ueta,ufai,xse,eta,fai,dx1,dx2,dx3,adiff,rr1;
+ double u1,u2,u3,sint[9];
+ double uxse,ueta,ufai,xse,eta,fai,adiff;
double prod1,unorm,twodiff;
- const int dims=E->mesh.nsd;
- const int dofs=E->mesh.dof;
- const int lev=E->mesh.levmax;
- const int nno=E->lmesh.nno;
- const int ends=enodes[E->mesh.nsd];
- const int vpts=vpoints[E->mesh.nsd];
-
ienm=E->ien[m][el].node;
twodiff = 2.0*diffusion;
@@ -390,64 +388,22 @@
uc3 += E->N.ppt[GNPINDEX(i,1)]*VV[3][i];
}
- dx1=0.25*(E->x[m][1][ienm[3]]+E->x[m][1][ienm[4]]
- +E->x[m][1][ienm[7]]+E->x[m][1][ienm[8]]
- -E->x[m][1][ienm[1]]-E->x[m][1][ienm[2]]
- -E->x[m][1][ienm[5]]-E->x[m][1][ienm[6]]);
- dx2=0.25*(E->x[m][2][ienm[3]]+E->x[m][2][ienm[4]]
- +E->x[m][2][ienm[7]]+E->x[m][2][ienm[8]]
- -E->x[m][2][ienm[1]]-E->x[m][2][ienm[2]]
- -E->x[m][2][ienm[5]]-E->x[m][2][ienm[6]]);
- dx3=0.25*(E->x[m][3][ienm[3]]+E->x[m][3][ienm[4]]
- +E->x[m][3][ienm[7]]+E->x[m][3][ienm[8]]
- -E->x[m][3][ienm[1]]-E->x[m][3][ienm[2]]
- -E->x[m][3][ienm[5]]-E->x[m][3][ienm[6]]);
- uxse = fabs(uc1*dx1+uc2*dx2+uc3*dx3);
+ uxse = fabs(uc1*E->eco[m][el].size[1]);
+ ueta = fabs(uc2*E->eco[m][el].size[2]);
+ ufai = fabs(uc3*E->eco[m][el].size[3]);
- dx1=0.25*(E->x[m][1][ienm[2]]+E->x[m][1][ienm[3]]
- +E->x[m][1][ienm[6]]+E->x[m][1][ienm[7]]
- -E->x[m][1][ienm[1]]-E->x[m][1][ienm[4]]
- -E->x[m][1][ienm[5]]-E->x[m][1][ienm[8]]);
- dx2=0.25*(E->x[m][2][ienm[2]]+E->x[m][2][ienm[3]]
- +E->x[m][2][ienm[6]]+E->x[m][2][ienm[7]]
- -E->x[m][2][ienm[1]]-E->x[m][2][ienm[4]]
- -E->x[m][2][ienm[5]]-E->x[m][2][ienm[8]]);
- dx3=0.25*(E->x[m][3][ienm[2]]+E->x[m][3][ienm[3]]
- +E->x[m][3][ienm[6]]+E->x[m][3][ienm[7]]
- -E->x[m][3][ienm[1]]-E->x[m][3][ienm[4]]
- -E->x[m][3][ienm[5]]-E->x[m][3][ienm[8]]);
- ueta = fabs(uc1*dx1+uc2*dx2+uc3*dx3);
-
- dx1=0.25*(E->x[m][1][ienm[5]]+E->x[m][1][ienm[6]]
- +E->x[m][1][ienm[7]]+E->x[m][1][ienm[8]]
- -E->x[m][1][ienm[1]]-E->x[m][1][ienm[2]]
- -E->x[m][1][ienm[3]]-E->x[m][1][ienm[4]]);
- dx2=0.25*(E->x[m][2][ienm[5]]+E->x[m][2][ienm[6]]
- +E->x[m][2][ienm[7]]+E->x[m][2][ienm[8]]
- -E->x[m][2][ienm[1]]-E->x[m][2][ienm[2]]
- -E->x[m][2][ienm[3]]-E->x[m][2][ienm[4]]);
- dx3=0.25*(E->x[m][3][ienm[5]]+E->x[m][3][ienm[6]]
- +E->x[m][3][ienm[7]]+E->x[m][3][ienm[8]]
- -E->x[m][3][ienm[1]]-E->x[m][3][ienm[2]]
- -E->x[m][3][ienm[3]]-E->x[m][3][ienm[4]]);
- ufai = fabs(uc1*dx1+uc2*dx2+uc3*dx3);
-
-/* xse = (uxse>twodiff)? (1.0-twodiff/uxse):0.0;
+ xse = (uxse>twodiff)? (1.0-twodiff/uxse):0.0;
eta = (ueta>twodiff)? (1.0-twodiff/ueta):0.0;
fai = (ufai>twodiff)? (1.0-twodiff/ufai):0.0;
-*/
- aa = 2.0*uxse/twodiff;
- bb = 2.0*ueta/twodiff;
- cc = 2.0*ufai/twodiff;
- xse = (1.0+exp(-aa))/(1.0-exp(-aa))-twodiff/uxse;
- eta = (1.0+exp(-bb))/(1.0-exp(-bb))-twodiff/ueta;
- fai = (1.0+exp(-cc))/(1.0-exp(-cc))-twodiff/ufai;
unorm = uc1*uc1 + uc2*uc2 + uc3*uc3;
adiff = (unorm>0.000001)?( (uxse*xse+ueta*eta+ufai*fai)/(2.0*unorm) ):0.0;
+ for(i=1;i<=VPOINTS3D;i++)
+ sint[i] = rtf[3][i]/sin(rtf[1][i]);
+
for(i=1;i<=VPOINTS3D;i++) {
u1 = u2 = u3 = 0.0;
for(j=1;j<=ENODES3D;j++) /* this line heavily used */ {
@@ -457,9 +413,9 @@
}
for(j=1;j<=ENODES3D;j++) {
- prod1 = (u1 * GNx->vpt[GNVXINDEX(0,j,i)] +
- u2 * GNx->vpt[GNVXINDEX(1,j,i)] +
- u3 * GNx->vpt[GNVXINDEX(2,j,i)] ) ;
+ prod1 = (u1 * GNx->vpt[GNVXINDEX(0,j,i)]*rtf[3][i] +
+ u2 * GNx->vpt[GNVXINDEX(1,j,i)]*sint[i] +
+ u3 * GNx->vpt[GNVXINDEX(2,j,i)] ) ;
PG->vpt[GNVINDEX(j,i)] = E->N.vpt[GNVINDEX(j,i)] + adiff * prod1;
}
@@ -475,7 +431,7 @@
Used to correct the Tdot term.
========================================= */
-void element_residual(E,el,PG,GNx,dOmega,VV,field,fielddot,Q0,Eres,diff,BC,FLAGS,m)
+void element_residual(E,el,PG,GNx,dOmega,VV,field,fielddot,Q0,Eres,rtf,diff,BC,FLAGS,m)
struct All_variables *E;
int el,m;
struct Shape_function PG;
@@ -485,6 +441,7 @@
double **field,**fielddot;
struct SOURCES Q0;
double Eres[9];
+ double rtf[4][9];
double diff;
float **BC;
unsigned int **FLAGS;
@@ -493,7 +450,7 @@
int i,j,a,k,node,nodes[5],d,aid,back_front,onedfns;
double Q;
double dT[9];
- double tx1[9],tx2[9],tx3[9];
+ double tx1[9],tx2[9],tx3[9],sint[9];
double v1[9],v2[9],v3[9];
double adv_dT,t2[4];
double T,DT;
@@ -520,6 +477,9 @@
v3[i] = tx3[i]= 0.0;
}
+ for(i=1;i<=vpts;i++)
+ sint[i] = rtf[3][i]/sin(rtf[1][i]);
+
for(j=1;j<=ends;j++) {
node = E->ien[m][el].node[j];
T = field[m][node];
@@ -530,8 +490,8 @@
for(i=1;i<=vpts;i++) {
dT[i] += DT * E->N.vpt[GNVINDEX(j,i)];
- tx1[i] += GNx.vpt[GNVXINDEX(0,j,i)] * T;
- tx2[i] += GNx.vpt[GNVXINDEX(1,j,i)] * T;
+ tx1[i] += GNx.vpt[GNVXINDEX(0,j,i)] * T * rtf[3][i];
+ tx2[i] += GNx.vpt[GNVXINDEX(1,j,i)] * T * sint[i];
tx3[i] += GNx.vpt[GNVXINDEX(2,j,i)] * T;
sfn = E->N.vpt[GNVINDEX(j,i)];
v1[i] += VV[1][j] * sfn;
@@ -557,8 +517,8 @@
Eres[j] -=
PG.vpt[GNVINDEX(j,i)] * dOmega.vpt[i]
* (dT[i] - Q + v1[i]*tx1[i] + v2[i]*tx2[i] + v3[i]*tx3[i])
- + diff*dOmega.vpt[i] * (GNx.vpt[GNVXINDEX(0,j,i)]*tx1[i] +
- GNx.vpt[GNVXINDEX(1,j,i)]*tx2[i] +
+ + diff*dOmega.vpt[i] * (GNx.vpt[GNVXINDEX(0,j,i)]*tx1[i]*rtf[3][i] +
+ GNx.vpt[GNVXINDEX(1,j,i)]*tx2[i]*sint[i] +
GNx.vpt[GNVXINDEX(2,j,i)]*tx3[i] );
}
}
More information about the cig-commits
mailing list