[cig-commits] r15398 - mc/3D/CitcomS/trunk/lib
leif at geodynamics.org
leif at geodynamics.org
Mon Jun 29 15:10:40 PDT 2009
Author: leif
Date: 2009-06-29 15:10:40 -0700 (Mon, 29 Jun 2009)
New Revision: 15398
Modified:
mc/3D/CitcomS/trunk/lib/multigrid_kernel.cu
Log:
For my own sanity, simplified code assuming the following:
E->sphere.caps_per_proc == 1
E->parallel.nproc == 1
E->mesh.levmax == 0
E->mesh.nsd == 3
Modified: mc/3D/CitcomS/trunk/lib/multigrid_kernel.cu
===================================================================
--- mc/3D/CitcomS/trunk/lib/multigrid_kernel.cu 2009-06-29 20:24:30 UTC (rev 15397)
+++ mc/3D/CitcomS/trunk/lib/multigrid_kernel.cu 2009-06-29 22:10:40 UTC (rev 15398)
@@ -5,21 +5,49 @@
#include "element_definitions.h"
+enum {
+ LEVEL = 0,
+ CAPS_PER_PROC = 1,
+ M = 1, /* cap # */
+ NSD = 3, /* Spatial extent: 3d */
+ MAX_EQN = NSD*14,
+};
+
+
+struct Some_variables {
+ int num_zero_resid;
+ int *zero_resid;
+
+ struct /*MESH_DATA*/ {
+ int NEQ;
+ int NNO;
+ } lmesh;
+
+ struct ID *ID;
+
+ higher_precision *Eqn_k1, *Eqn_k2, *Eqn_k3;
+ int *Node_map;
+
+ double *BI;
+
+ double *temp, *temp1;
+ unsigned int *NODE;
+};
+
+
/*------------------------------------------------------------------------*/
/* from BC_util.c */
__device__ void strip_bcs_from_residual(
- struct All_variables *E,
- double **Res,
- int level
+ struct Some_variables *E,
+ double *Res
)
{
- int m,i;
+ int i;
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- if (E->num_zero_resid[level][m])
- for(i=1;i<=E->num_zero_resid[level][m];i++)
- Res[m][E->zero_resid[level][m][i]] = 0.0;
+ if (E->num_zero_resid)
+ for(i=1;i<=E->num_zero_resid;i++)
+ Res[E->zero_resid[i]] = 0.0;
return;
}
@@ -29,13 +57,12 @@
/* from Element_calculations.c */
__device__ void n_assemble_del2_u(
- struct All_variables *E,
- double **u, double **Au,
- int level,
+ struct Some_variables *E,
+ double *u, double *Au,
int strip_bcs
)
{
- int m, e,i;
+ int e,i;
int eqn1,eqn2,eqn3;
double UU,U1,U2,U3;
@@ -43,50 +70,43 @@
int *C;
higher_precision *B1,*B2,*B3;
- const int neq=E->lmesh.NEQ[level];
- const int nno=E->lmesh.NNO[level];
- const int dims=E->mesh.nsd;
- const int max_eqn = dims*14;
+ const int neq=E->lmesh.NEQ;
+ const int nno=E->lmesh.NNO;
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ for(e=0;e<=neq;e++)
+ Au[e]=0.0;
- for(e=0;e<=neq;e++)
- Au[m][e]=0.0;
+ u[neq] = 0.0;
- u[m][neq] = 0.0;
+ for(e=1;e<=nno;e++) {
- for(e=1;e<=nno;e++) {
+ eqn1=E->ID[e].doff[1];
+ eqn2=E->ID[e].doff[2];
+ eqn3=E->ID[e].doff[3];
- eqn1=E->ID[level][m][e].doff[1];
- eqn2=E->ID[level][m][e].doff[2];
- eqn3=E->ID[level][m][e].doff[3];
+ U1 = u[eqn1];
+ U2 = u[eqn2];
+ U3 = u[eqn3];
- U1 = u[m][eqn1];
- U2 = u[m][eqn2];
- U3 = u[m][eqn3];
+ C=E->Node_map + (e-1)*MAX_EQN;
+ B1=E->Eqn_k1+(e-1)*MAX_EQN;
+ B2=E->Eqn_k2+(e-1)*MAX_EQN;
+ B3=E->Eqn_k3+(e-1)*MAX_EQN;
- C=E->Node_map[level][m] + (e-1)*max_eqn;
- B1=E->Eqn_k1[level][m]+(e-1)*max_eqn;
- B2=E->Eqn_k2[level][m]+(e-1)*max_eqn;
- B3=E->Eqn_k3[level][m]+(e-1)*max_eqn;
+ for(i=3;i<MAX_EQN;i++) {
+ UU = u[C[i]];
+ Au[eqn1] += B1[i]*UU;
+ Au[eqn2] += B2[i]*UU;
+ Au[eqn3] += B3[i]*UU;
+ }
+ for(i=0;i<MAX_EQN;i++)
+ Au[C[i]] += B1[i]*U1+B2[i]*U2+B3[i]*U3;
- for(i=3;i<max_eqn;i++) {
- UU = u[m][C[i]];
- Au[m][eqn1] += B1[i]*UU;
- Au[m][eqn2] += B2[i]*UU;
- Au[m][eqn3] += B3[i]*UU;
- }
- for(i=0;i<max_eqn;i++)
- Au[m][C[i]] += B1[i]*U1+B2[i]*U2+B3[i]*U3;
+ } /* end for e */
- } /* end for e */
- } /* end for m */
-
- (E->solver.exchange_id_d)(E, Au, level);
-
if (strip_bcs)
- strip_bcs_from_residual(E,Au,level);
+ strip_bcs_from_residual(E,Au);
return;
}
@@ -96,17 +116,16 @@
/* based on the function from General_matrix_functions.c */
__global__ void do_gauss_seidel(
- struct All_variables *E,
- double **d0,
- double **F, double **Ad,
+ struct Some_variables *E,
+ double *d0,
+ double *F, double *Ad,
double acc,
int *cycles,
- int level,
int guess
)
{
- int count,i,j,m,steps;
+ int count,i,j,steps;
int *C;
int eqn1,eqn2,eqn3;
@@ -115,121 +134,104 @@
higher_precision *B1,*B2,*B3;
- const int dims=E->mesh.nsd;
- const int neq=E->lmesh.NEQ[level];
- const int max_eqn=14*dims;
+ const int neq=E->lmesh.NEQ;
const double zeroo = 0.0;
steps=*cycles;
if(guess) {
- n_assemble_del2_u(E,d0,Ad,level,1);
+ n_assemble_del2_u(E,d0,Ad,1);
}
else
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++) {
- d0[m][i]=Ad[m][i]=zeroo;
- }
+ for(i=0;i<neq;i++) {
+ d0[i]=Ad[i]=zeroo;
+ }
count = 0;
while (count < steps) {
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(j=0;j<=E->lmesh.NEQ[level];j++)
- E->temp[m][j] = zeroo;
+ for(j=0;j<=E->lmesh.NEQ;j++)
+ E->temp[j] = zeroo;
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- Ad[m][neq] = zeroo;
+ Ad[neq] = zeroo;
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.NNO[level];i++)
- if(E->NODE[level][m][i] & OFFSIDE) {
+ for(i=1;i<=E->lmesh.NNO;i++)
+ if(E->NODE[i] & OFFSIDE) {
- eqn1=E->ID[level][m][i].doff[1];
- eqn2=E->ID[level][m][i].doff[2];
- eqn3=E->ID[level][m][i].doff[3];
- E->temp[m][eqn1] = (F[m][eqn1] - Ad[m][eqn1])*E->BI[level][m][eqn1];
- E->temp[m][eqn2] = (F[m][eqn2] - Ad[m][eqn2])*E->BI[level][m][eqn2];
- E->temp[m][eqn3] = (F[m][eqn3] - Ad[m][eqn3])*E->BI[level][m][eqn3];
- E->temp1[m][eqn1] = Ad[m][eqn1];
- E->temp1[m][eqn2] = Ad[m][eqn2];
- E->temp1[m][eqn3] = Ad[m][eqn3];
- }
+ eqn1=E->ID[i].doff[1];
+ eqn2=E->ID[i].doff[2];
+ eqn3=E->ID[i].doff[3];
+ E->temp[eqn1] = (F[eqn1] - Ad[eqn1])*E->BI[eqn1];
+ E->temp[eqn2] = (F[eqn2] - Ad[eqn2])*E->BI[eqn2];
+ E->temp[eqn3] = (F[eqn3] - Ad[eqn3])*E->BI[eqn3];
+ E->temp1[eqn1] = Ad[eqn1];
+ E->temp1[eqn2] = Ad[eqn2];
+ E->temp1[eqn3] = Ad[eqn3];
+ }
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.NNO[level];i++) {
+ for(i=1;i<=E->lmesh.NNO;i++) {
- eqn1=E->ID[level][m][i].doff[1];
- eqn2=E->ID[level][m][i].doff[2];
- eqn3=E->ID[level][m][i].doff[3];
- C=E->Node_map[level][m]+(i-1)*max_eqn;
- B1=E->Eqn_k1[level][m]+(i-1)*max_eqn;
- B2=E->Eqn_k2[level][m]+(i-1)*max_eqn;
- B3=E->Eqn_k3[level][m]+(i-1)*max_eqn;
+ eqn1=E->ID[i].doff[1];
+ eqn2=E->ID[i].doff[2];
+ eqn3=E->ID[i].doff[3];
+ C=E->Node_map+(i-1)*MAX_EQN;
+ B1=E->Eqn_k1+(i-1)*MAX_EQN;
+ B2=E->Eqn_k2+(i-1)*MAX_EQN;
+ B3=E->Eqn_k3+(i-1)*MAX_EQN;
- /* Ad on boundaries differs after the following operation, but
- no communications are needed yet, because boundary Ad will
- not be used for the G-S iterations for interior nodes */
+ /* Ad on boundaries differs after the following operation, but
+ no communications are needed yet, because boundary Ad will
+ not be used for the G-S iterations for interior nodes */
- for(j=3;j<max_eqn;j++) {
- UU = E->temp[m][C[j]];
- Ad[m][eqn1] += B1[j]*UU;
- Ad[m][eqn2] += B2[j]*UU;
- Ad[m][eqn3] += B3[j]*UU;
- }
+ for(j=3;j<MAX_EQN;j++) {
+ UU = E->temp[C[j]];
+ Ad[eqn1] += B1[j]*UU;
+ Ad[eqn2] += B2[j]*UU;
+ Ad[eqn3] += B3[j]*UU;
+ }
- if (!(E->NODE[level][m][i]&OFFSIDE)) {
- E->temp[m][eqn1] = (F[m][eqn1] - Ad[m][eqn1])*E->BI[level][m][eqn1];
- E->temp[m][eqn2] = (F[m][eqn2] - Ad[m][eqn2])*E->BI[level][m][eqn2];
- E->temp[m][eqn3] = (F[m][eqn3] - Ad[m][eqn3])*E->BI[level][m][eqn3];
- }
+ if (!(E->NODE[i]&OFFSIDE)) {
+ E->temp[eqn1] = (F[eqn1] - Ad[eqn1])*E->BI[eqn1];
+ E->temp[eqn2] = (F[eqn2] - Ad[eqn2])*E->BI[eqn2];
+ E->temp[eqn3] = (F[eqn3] - Ad[eqn3])*E->BI[eqn3];
+ }
- /* Ad on boundaries differs after the following operation */
- for(j=0;j<max_eqn;j++)
- Ad[m][C[j]] += B1[j]*E->temp[m][eqn1]
- + B2[j]*E->temp[m][eqn2]
- + B3[j]*E->temp[m][eqn3];
+ /* Ad on boundaries differs after the following operation */
+ for(j=0;j<MAX_EQN;j++)
+ Ad[C[j]] += B1[j]*E->temp[eqn1]
+ + B2[j]*E->temp[eqn2]
+ + B3[j]*E->temp[eqn3];
- d0[m][eqn1] += E->temp[m][eqn1];
- d0[m][eqn2] += E->temp[m][eqn2];
- d0[m][eqn3] += E->temp[m][eqn3];
- }
+ d0[eqn1] += E->temp[eqn1];
+ d0[eqn2] += E->temp[eqn2];
+ d0[eqn3] += E->temp[eqn3];
+ }
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.NNO[level];i++)
- if(E->NODE[level][m][i] & OFFSIDE) {
- eqn1=E->ID[level][m][i].doff[1];
- eqn2=E->ID[level][m][i].doff[2];
- eqn3=E->ID[level][m][i].doff[3];
- Ad[m][eqn1] -= E->temp1[m][eqn1];
- Ad[m][eqn2] -= E->temp1[m][eqn2];
- Ad[m][eqn3] -= E->temp1[m][eqn3];
- }
+ for(i=1;i<=E->lmesh.NNO;i++)
+ if(E->NODE[i] & OFFSIDE) {
+ eqn1=E->ID[i].doff[1];
+ eqn2=E->ID[i].doff[2];
+ eqn3=E->ID[i].doff[3];
+ Ad[eqn1] -= E->temp1[eqn1];
+ Ad[eqn2] -= E->temp1[eqn2];
+ Ad[eqn3] -= E->temp1[eqn3];
+ }
- (E->solver.exchange_id_d)(E, Ad, level);
+ for(i=1;i<=E->lmesh.NNO;i++)
+ if(E->NODE[i] & OFFSIDE) {
+ eqn1=E->ID[i].doff[1];
+ eqn2=E->ID[i].doff[2];
+ eqn3=E->ID[i].doff[3];
+ Ad[eqn1] += E->temp1[eqn1];
+ Ad[eqn2] += E->temp1[eqn2];
+ Ad[eqn3] += E->temp1[eqn3];
+ }
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.NNO[level];i++)
- if(E->NODE[level][m][i] & OFFSIDE) {
- eqn1=E->ID[level][m][i].doff[1];
- eqn2=E->ID[level][m][i].doff[2];
- eqn3=E->ID[level][m][i].doff[3];
- Ad[m][eqn1] += E->temp1[m][eqn1];
- Ad[m][eqn2] += E->temp1[m][eqn2];
- Ad[m][eqn3] += E->temp1[m][eqn3];
- }
-
count++;
- /* for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++) {
- F[m][i] -= Ad[m][i];
- Ad[m][i] = 0.0;
- }
- */
}
*cycles=count;
@@ -240,6 +242,19 @@
/*------------------------------------------------------------------------*/
+static void assert_assumptions(struct All_variables *E, int level) {
+
+ assert(E->control.NMULTIGRID);
+
+ assert(E->sphere.caps_per_proc == CAPS_PER_PROC);
+
+ assert(E->mesh.nsd == NSD);
+ assert(E->mesh.levmax == LEVEL);
+ assert(level == LEVEL);
+
+ assert(E->parallel.nproc == 1);
+}
+
extern "C" void gauss_seidel(
struct All_variables *E,
double **d0,
@@ -250,5 +265,31 @@
int guess
)
{
+ struct Some_variables kE;
+
+ assert_assumptions(E, level);
+
+ /* initialize 'Some_variables' with 'All_variables' */
+
+ kE.num_zero_resid = E->num_zero_resid[LEVEL][M];
+ kE.zero_resid = E->zero_resid[LEVEL][M];
+
+ kE.lmesh.NEQ = E->lmesh.NEQ[LEVEL];
+ kE.lmesh.NNO = E->lmesh.NNO[LEVEL];
+
+ kE.ID = E->ID[LEVEL][M];
+
+ kE.Eqn_k1 = E->Eqn_k1[LEVEL][M];
+ kE.Eqn_k2 = E->Eqn_k2[LEVEL][M];
+ kE.Eqn_k3 = E->Eqn_k3[LEVEL][M];
+ kE.Node_map = E->Node_map[LEVEL][M];
+
+ kE.BI = E->BI[LEVEL][M];
+
+ kE.temp = E->temp[M];
+ kE.temp1 = E->temp1[M];
+
+ kE.NODE = E->NODE[LEVEL][M];
+
/* XXX */
}
More information about the CIG-COMMITS
mailing list