[cig-commits] r14809 - mc/3D/CitcomS/trunk/lib
leif at geodynamics.org
leif at geodynamics.org
Tue Apr 28 15:53:39 PDT 2009
Author: leif
Date: 2009-04-28 15:53:38 -0700 (Tue, 28 Apr 2009)
New Revision: 14809
Modified:
mc/3D/CitcomS/trunk/lib/cgrad_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/cgrad_kernel.cu
===================================================================
--- mc/3D/CitcomS/trunk/lib/cgrad_kernel.cu 2009-04-28 21:28:34 UTC (rev 14808)
+++ mc/3D/CitcomS/trunk/lib/cgrad_kernel.cu 2009-04-28 22:53:38 UTC (rev 14809)
@@ -5,104 +5,80 @@
#include "global_defs.h"
-#ifndef __CUDACC__
-#define __device__
-#define __global__
-#define __host__
-#else
-/* XXX */
+#ifdef __CUDACC__
#undef assert
#define assert(c)
+#else
+#define __device__ static inline
+#define __global__ static
+#define __host__
#endif
-/*------------------------------------------------------------------------*/
-/* from element_definitions.h */
+enum {
+ LEVEL = 0,
+ CAPS_PER_PROC = 1,
+ M = 1, /* cap # */
+ NSD = 3, /* Spatial extent: 3d */
+ MAX_EQN = NSD*14,
+ ENODES = 8, /* enodes[NSD] */
+ LOC_MAT_SIZE = 24, /* loc_mat_size[NSD] */
+};
-__device__ static const int enodes[] = {0,2,4,8};
-__device__ static const int loc_mat_size[] = {0,4,8,24};
-
-/*------------------------------------------------------------------------*/
-/* from Regional_parallel_related.c */
-
-__device__ void regional_exchange_id_d(
- struct All_variables *E,
- double **U,
- int lev,
- double *memory
- )
-{
-
- int j,m,k;
- double *S[27],*R[27];
- int dimofk;
+struct Some_variables {
+ int num_zero_resid;
+ int *zero_resid;
-#if 0 /* XXX */
- MPI_Status status;
-#endif
+ struct /*MESH_DATA*/ {
+ int NEQ;
+ int NNO;
+ int NEL;
+ } lmesh;
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- for (k=1;k<=E->parallel.TNUM_PASS[lev][m];k++) {
- dimofk = 1+E->parallel.NUM_NEQ[lev][m].pass[k];
- S[k] = memory; memory += dimofk;
- R[k] = memory; memory += dimofk;
- }
- }
-
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- for (k=1;k<=E->parallel.TNUM_PASS[lev][m];k++) {
-
- for (j=1;j<=E->parallel.NUM_NEQ[lev][m].pass[k];j++)
- S[k][j-1] = U[m][ E->parallel.EXCHANGE_ID[lev][m][j].pass[k] ];
-
-#if 0 /* XXX */
- MPI_Sendrecv(S[k],E->parallel.NUM_NEQ[lev][m].pass[k],MPI_DOUBLE,
- E->parallel.PROCESSOR[lev][m].pass[k],1,
- R[k],E->parallel.NUM_NEQ[lev][m].pass[k],MPI_DOUBLE,
- E->parallel.PROCESSOR[lev][m].pass[k],1,
- E->parallel.world,&status);
-#endif
-
- for (j=1;j<=E->parallel.NUM_NEQ[lev][m].pass[k];j++)
- U[m][ E->parallel.EXCHANGE_ID[lev][m][j].pass[k] ] += R[k][j-1];
-
- } /* for k */
- } /* for m */ /* finish sending */
+ struct IEN *IEN;
+ struct ID *ID;
+ struct EK *elt_k;
- return;
-}
+ higher_precision *Eqn_k1, *Eqn_k2, *Eqn_k3;
+ int *Node_map;
+
+ double *BI;
+
+ struct /*CONTROL*/ {
+ int NASSEMBLE;
+
+ int v_steps_low;
+ int total_iteration_cycles;
+ int total_v_solver_calls;
+ } control;
+
+ /* temporary malloc'd memory */
+ double *memory;
+ int memoryDim;
+
+ /* outputs */
+
+ struct /*MONITOR*/ {
+ double momentum_residual;
+ } monitor;
+ int valid;
+};
-/*------------------------------------------------------------------------*/
-/* from Full_parallel_related.c */
-__device__ void full_exchange_id_d(
- struct All_variables *E,
- double **U,
- int lev,
- double *memory
- )
-{
- /* XXX */
-}
-
-
/*------------------------------------------------------------------------*/
/* 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;
+ for(i=1;i<=E->num_zero_resid;i++)
+ Res[E->zero_resid[i]] = 0.0;
return;
}
@@ -112,87 +88,72 @@
/* from Element_calculations.c */
__device__ void e_assemble_del2_u(
- struct All_variables *E,
- double **u, double **Au,
- int level,
- int strip_bcs,
- double *memory
+ struct Some_variables *E,
+ double *u, double *Au,
+ int strip_bcs
)
{
- int e,i,a,b,a1,a2,a3,ii,m,nodeb;
+ int e,i,a,b,a1,a2,a3,ii,nodeb;
- const int n=loc_mat_size[E->mesh.nsd];
- const int ends=enodes[E->mesh.nsd];
- const int dims=E->mesh.nsd;
- const int nel=E->lmesh.NEL[level];
- const int neq=E->lmesh.NEQ[level];
+ const int nel=E->lmesh.NEL;
+ const int neq=E->lmesh.NEQ;
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- for(i=0;i<neq;i++)
- Au[m][i] = 0.0;
+ for(i=0;i<neq;i++)
+ Au[i] = 0.0;
- for(e=1;e<=nel;e++) {
- for(a=1;a<=ends;a++) {
- ii = E->IEN[level][m][e].node[a];
- a1 = E->ID[level][m][ii].doff[1];
- a2 = E->ID[level][m][ii].doff[2];
- a3 = E->ID[level][m][ii].doff[3];
- for(b=1;b<=ends;b++) {
- nodeb = E->IEN[level][m][e].node[b];
- ii = (a*n+b)*dims-(dims*n+dims);
- /* i=1, j=1,2 */
- /* i=1, j=1,2,3 */
- Au[m][a1] +=
- E->elt_k[level][m][e].k[ii] *
- u[m][E->ID[level][m][nodeb].doff[1]]
- + E->elt_k[level][m][e].k[ii+1] *
- u[m][E->ID[level][m][nodeb].doff[2]]
- + E->elt_k[level][m][e].k[ii+2] *
- u[m][E->ID[level][m][nodeb].doff[3]];
- /* i=2, j=1,2,3 */
- Au[m][a2] +=
- E->elt_k[level][m][e].k[ii+n] *
- u[m][E->ID[level][m][nodeb].doff[1]]
- + E->elt_k[level][m][e].k[ii+n+1] *
- u[m][E->ID[level][m][nodeb].doff[2]]
- + E->elt_k[level][m][e].k[ii+n+2] *
- u[m][E->ID[level][m][nodeb].doff[3]];
- /* i=3, j=1,2,3 */
- Au[m][a3] +=
- E->elt_k[level][m][e].k[ii+n+n] *
- u[m][E->ID[level][m][nodeb].doff[1]]
- + E->elt_k[level][m][e].k[ii+n+n+1] *
- u[m][E->ID[level][m][nodeb].doff[2]]
- + E->elt_k[level][m][e].k[ii+n+n+2] *
- u[m][E->ID[level][m][nodeb].doff[3]];
+ for(e=1;e<=nel;e++) {
+ for(a=1;a<=ENODES;a++) {
+ ii = E->IEN[e].node[a];
+ a1 = E->ID[ii].doff[1];
+ a2 = E->ID[ii].doff[2];
+ a3 = E->ID[ii].doff[3];
+ for(b=1;b<=ENODES;b++) {
+ nodeb = E->IEN[e].node[b];
+ ii = (a*LOC_MAT_SIZE+b)*NSD-(NSD*LOC_MAT_SIZE+NSD);
+ /* i=1, j=1,2 */
+ /* i=1, j=1,2,3 */
+ Au[a1] +=
+ E->elt_k[e].k[ii] *
+ u[E->ID[nodeb].doff[1]]
+ + E->elt_k[e].k[ii+1] *
+ u[E->ID[nodeb].doff[2]]
+ + E->elt_k[e].k[ii+2] *
+ u[E->ID[nodeb].doff[3]];
+ /* i=2, j=1,2,3 */
+ Au[a2] +=
+ E->elt_k[e].k[ii+LOC_MAT_SIZE] *
+ u[E->ID[nodeb].doff[1]]
+ + E->elt_k[e].k[ii+LOC_MAT_SIZE+1] *
+ u[E->ID[nodeb].doff[2]]
+ + E->elt_k[e].k[ii+LOC_MAT_SIZE+2] *
+ u[E->ID[nodeb].doff[3]];
+ /* i=3, j=1,2,3 */
+ Au[a3] +=
+ E->elt_k[e].k[ii+LOC_MAT_SIZE+LOC_MAT_SIZE] *
+ u[E->ID[nodeb].doff[1]]
+ + E->elt_k[e].k[ii+LOC_MAT_SIZE+LOC_MAT_SIZE+1] *
+ u[E->ID[nodeb].doff[2]]
+ + E->elt_k[e].k[ii+LOC_MAT_SIZE+LOC_MAT_SIZE+2] *
+ u[E->ID[nodeb].doff[3]];
- } /* end for loop b */
- } /* end for loop a */
+ } /* end for loop b */
+ } /* end for loop a */
- } /* end for e */
- } /* end for m */
+ } /* end for e */
- if (0) { /* XXX */
- full_exchange_id_d(E, Au, level, memory);
- } else {
- regional_exchange_id_d(E, Au, level, memory);
- }
-
if(strip_bcs)
- strip_bcs_from_residual(E,Au,level);
+ strip_bcs_from_residual(E,Au);
return;
}
__device__ void n_assemble_del2_u(
- struct All_variables *E,
- double **u, double **Au,
- int level,
- int strip_bcs,
- double *memory
+ 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;
@@ -200,70 +161,61 @@
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;
+
+ /*
+ * Au = E->Eqn_k? * u
+ * where E->Eqn_k? is the sparse stiffness matrix
+ */
+ for(e=0;e<=neq;e++)
+ Au[e]=0.0;
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ u[neq] = 0.0;
- for(e=0;e<=neq;e++)
- Au[m][e]=0.0;
+ for(e=1;e<=nno;e++) {
- u[m][neq] = 0.0;
+ eqn1=E->ID[e].doff[1];
+ eqn2=E->ID[e].doff[2];
+ eqn3=E->ID[e].doff[3];
- for(e=1;e<=nno;e++) {
+ U1 = u[eqn1];
+ U2 = u[eqn2];
+ U3 = u[eqn3];
- eqn1=E->ID[level][m][e].doff[1];
- eqn2=E->ID[level][m][e].doff[2];
- eqn3=E->ID[level][m][e].doff[3];
+ 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;
- U1 = u[m][eqn1];
- U2 = u[m][eqn2];
- U3 = u[m][eqn3];
+ 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;
- 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[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 m */
+ } /* end for e */
- if (0) { /* XXX */
- full_exchange_id_d(E, Au, level, memory);
- } else {
- regional_exchange_id_d(E, Au, level, memory);
- }
-
if (strip_bcs)
- strip_bcs_from_residual(E,Au,level);
+ strip_bcs_from_residual(E,Au);
return;
}
__device__ void assemble_del2_u(
- struct All_variables *E,
- double **u, double **Au,
- int level,
- int strip_bcs,
- double *memory
+ struct Some_variables *E,
+ double *u, double *Au,
+ int strip_bcs
)
{
- if(E->control.NMULTIGRID||E->control.NASSEMBLE)
- n_assemble_del2_u(E,u,Au,level,strip_bcs, memory);
+ if (E->control.NASSEMBLE)
+ n_assemble_del2_u(E,u,Au,strip_bcs);
else
- e_assemble_del2_u(E,u,Au,level,strip_bcs, memory);
+ e_assemble_del2_u(E,u,Au,strip_bcs);
return;
}
@@ -273,35 +225,20 @@
/* from Global_operations.c */
__device__ double global_vdot(
- struct All_variables *E,
- double **A, double **B,
- int lev
+ struct Some_variables *E,
+ double *A, double *B
)
{
- int m,i,neq;
- double prod, temp,temp1;
+ int i,neq;
+ double prod;
- temp = 0.0;
prod = 0.0;
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- neq=E->lmesh.NEQ[lev];
- temp1 = 0.0;
- for (i=0;i<neq;i++)
- temp += A[m][i]*B[m][i];
+ neq=E->lmesh.NEQ;
+ for (i=0;i<neq;i++)
+ prod += A[i]*B[i];
- for (i=1;i<=E->parallel.Skip_neq[lev][m];i++)
- temp1 += A[m][E->parallel.Skip_id[lev][m][i]]*B[m][E->parallel.Skip_id[lev][m][i]];
-
- temp -= temp1;
-
- }
-
-#if 0 /* XXX */
- MPI_Allreduce(&temp, &prod,1,MPI_DOUBLE,MPI_SUM,E->parallel.world);
-#endif
-
- return (prod);
+ return prod;
}
@@ -309,96 +246,87 @@
/* from General_matix_functions.c */
__device__ double conj_grad(
- struct All_variables *E,
- double **d0,
- double **F,
+ struct Some_variables *E,
+ double *d0,
+ double *F,
double acc,
- int *cycles,
- int level,
- double *memory
+ int *cycles
)
{
- double *r0[NCS],*r1[NCS],*r2[NCS];
- double *z0[NCS],*z1[NCS];
- double *p1[NCS],*p2[NCS];
- double *Ap[NCS];
- double *shuffle[NCS];
+ double *r0,*r1,*r2;
+ double *z0,*z1;
+ double *p1,*p2;
+ double *Ap;
+ double *shuffle;
- int m,count,i,steps;
+ int count,i,steps;
double residual;
double alpha,beta,dotprod,dotr1z1,dotr0z0;
+
+ double *memory = E->memory;
- const int mem_lev=E->mesh.levmax;
- const int high_neq = E->lmesh.NEQ[level];
+ const int neq = E->lmesh.NEQ;
steps = *cycles;
- for(m=1;m<=E->sphere.caps_per_proc;m++) {
- r0[m] = memory; memory += E->lmesh.NEQ[mem_lev];
- r1[m] = memory; memory += E->lmesh.NEQ[mem_lev];
- r2[m] = memory; memory += E->lmesh.NEQ[mem_lev];
- z0[m] = memory; memory += E->lmesh.NEQ[mem_lev];
- z1[m] = memory; memory += E->lmesh.NEQ[mem_lev];
- p1[m] = memory; memory += (1+E->lmesh.NEQ[mem_lev]);
- p2[m] = memory; memory += (1+E->lmesh.NEQ[mem_lev]);
- Ap[m] = memory; memory += (1+E->lmesh.NEQ[mem_lev]);
+ r0 = memory; memory += neq;
+ r1 = memory; memory += neq;
+ r2 = memory; memory += neq;
+ z0 = memory; memory += neq;
+ z1 = memory; memory += neq;
+ p1 = memory; memory += (1+neq);
+ p2 = memory; memory += (1+neq);
+ Ap = memory; memory += (1+neq);
+ assert(memory == E->memory + E->memoryDim);
+
+ for(i=0;i<neq;i++) {
+ r1[i] = F[i];
+ d0[i] = 0.0;
}
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<high_neq;i++) {
- r1[m][i] = F[m][i];
- d0[m][i] = 0.0;
- }
+ residual = sqrt(global_vdot(E,r1,r1));
- residual = sqrt(global_vdot(E,r1,r1,level));
-
assert(residual != 0.0 /* initial residual for CG = 0.0 */);
count = 0;
while (((residual > acc) && (count < steps)) || count == 0) {
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<high_neq;i++)
- z1[m][i] = E->BI[level][m][i] * r1[m][i];
+ for(i=0;i<neq;i++)
+ z1[i] = E->BI[i] * r1[i];
- dotr1z1 = global_vdot(E,r1,z1,level);
+ dotr1z1 = global_vdot(E,r1,z1);
if (0==count)
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<high_neq;i++)
- p2[m][i] = z1[m][i];
+ for(i=0;i<neq;i++)
+ p2[i] = z1[i];
else {
assert(dotr0z0 != 0.0 /* in head of conj_grad */);
beta = dotr1z1/dotr0z0;
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<high_neq;i++)
- p2[m][i] = z1[m][i] + beta * p1[m][i];
+ for(i=0;i<neq;i++)
+ p2[i] = z1[i] + beta * p1[i];
}
dotr0z0 = dotr1z1;
- assemble_del2_u(E,p2,Ap,level,1,memory);
+ assemble_del2_u(E,p2,Ap,1);
- dotprod=global_vdot(E,p2,Ap,level);
+ dotprod=global_vdot(E,p2,Ap);
if(0.0==dotprod)
alpha=1.0e-3;
else
alpha = dotr1z1/dotprod;
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<high_neq;i++) {
- d0[m][i] += alpha * p2[m][i];
- r2[m][i] = r1[m][i] - alpha * Ap[m][i];
- }
+ for(i=0;i<neq;i++) {
+ d0[i] += alpha * p2[i];
+ r2[i] = r1[i] - alpha * Ap[i];
+ }
- residual = sqrt(global_vdot(E,r2,r2,level));
+ residual = sqrt(global_vdot(E,r2,r2));
- for(m=1;m<=E->sphere.caps_per_proc;m++) {
- shuffle[m] = r0[m]; r0[m] = r1[m]; r1[m] = r2[m]; r2[m] = shuffle[m];
- shuffle[m] = z0[m]; z0[m] = z1[m]; z1[m] = shuffle[m];
- shuffle[m] = p1[m]; p1[m] = p2[m]; p2[m] = shuffle[m];
- }
+ shuffle = r0; r0 = r1; r1 = r2; r2 = shuffle;
+ shuffle = z0; z0 = z1; z1 = shuffle;
+ shuffle = p1; p1 = p2; p2 = shuffle;
count++;
/* end of while-loop */
@@ -407,44 +335,37 @@
*cycles=count;
- strip_bcs_from_residual(E,d0,level);
+ strip_bcs_from_residual(E,d0);
- return(residual);
+ return residual;
}
__global__ void solve_del2_u(
- struct All_variables *E,
- double **d0,
- double **F,
- double acc,
- int high_lev,
- double *memory,
- int *valid
+ struct Some_variables *E,
+ double *d0,
+ double *F,
+ double acc
)
{
int count,cycles;
- int i, neq, m;
+ int i, neq;
double residual;
- neq = E->lmesh.NEQ[high_lev];
+ neq = E->lmesh.NEQ;
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++) {
- d0[m][i] = 0.0;
- }
+ for(i=0;i<neq;i++) {
+ d0[i] = 0.0;
+ }
- residual = sqrt(global_vdot(E,F,F,high_lev));
+ residual = sqrt(global_vdot(E,F,F));
count = 0;
- assert(!E->control.NMULTIGRID);
- /* conjugate gradient solution */
-
cycles = E->control.v_steps_low;
- residual = conj_grad(E,d0,F,acc,&cycles,high_lev,memory);
- *valid = (residual < acc)? 1:0;
+ residual = conj_grad(E,d0,F,acc,&cycles);
+ E->valid = (residual < acc)? 1:0;
count++;
@@ -458,36 +379,78 @@
/*------------------------------------------------------------------------*/
-__host__ int main() {
-#ifndef __CUDACC__
- if (0) {
- struct All_variables *E;
- const int levmax = E->mesh.levmax;
- int m,k;
-
- double *memory;
- int memoryDim;
- int valid;
+static void assert_assumptions(struct All_variables *E, int high_lev) {
+
+ assert(!E->control.NMULTIGRID);
+
+ assert(E->sphere.caps_per_proc == CAPS_PER_PROC);
+
+ assert(E->mesh.nsd == NSD);
+ assert(E->mesh.levmax == LEVEL);
+ assert(high_lev == LEVEL);
+
+ assert(E->parallel.nproc == 1);
+ assert(E->parallel.TNUM_PASS[LEVEL][M] == 0);
+ assert(E->parallel.Skip_neq[LEVEL][M] == 0);
+}
- memoryDim = 0;
- /* conj_grad */
- memoryDim += E->sphere.caps_per_proc *
- (5 * E->lmesh.NEQ[levmax] + /* r0,r1,r2,z0,z1 */
- 3 * (1+E->lmesh.NEQ[levmax]) /* p1,p2,Ap */
- );
- /* regional_exchange_id_d */
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- for (k=1;k<=E->parallel.TNUM_PASS[levmax][m];k++) {
- memoryDim += 2 * (1+E->parallel.NUM_NEQ[levmax][m].pass[k]); /* S,R */
- }
- }
+__host__ int launch_solve_del2_u(
+ struct All_variables *E,
+ double **d0,
+ double **F,
+ double acc,
+ int high_lev
+ )
+{
+ struct Some_variables kE;
+
+ assert_assumptions(E, high_lev);
+
+ /* initialize 'Some_variables' with 'All_variables' */
+
+ kE.num_zero_resid = E->num_zero_resid[LEVEL][M];
+ kE.zero_resid = E->zero_resid[LEVEL][M];
- memory = (double *)malloc(memoryDim*sizeof(double));
+ kE.lmesh.NEQ = E->lmesh.NEQ[LEVEL];
+ kE.lmesh.NNO = E->lmesh.NNO[LEVEL];
+ kE.lmesh.NEL = E->lmesh.NEL[LEVEL];
- solve_del2_u(0,0,0,0,0,memory,&valid);
+ kE.IEN = E->IEN[LEVEL][M];
+ kE.ID = E->ID[LEVEL][M];
+ kE.elt_k = E->elt_k[LEVEL][M];
- free(memory);
- }
-#endif
- return 0;
+ 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.control.NASSEMBLE = E->control.NASSEMBLE;
+ kE.control.v_steps_low = E->control.v_steps_low;
+ kE.control.total_iteration_cycles = E->control.total_iteration_cycles; /* in/out */
+ kE.control.total_v_solver_calls = E->control.total_v_solver_calls; /* in/out */
+
+ /* allocate temporary memory */
+ kE.memoryDim = 0;
+ kE.memoryDim += 5 * E->lmesh.NEQ[LEVEL] + /* r0,r1,r2,z0,z1 */
+ 3 * (1+E->lmesh.NEQ[LEVEL]) /* p1,p2,Ap */
+ ;
+ kE.memory = (double *)malloc(kE.memoryDim*sizeof(double));
+
+ /* zero outputs */
+ kE.monitor.momentum_residual = 0.0;
+ kE.valid = 0;
+
+ /* solve */
+ solve_del2_u(&kE, d0[M], F[M], acc);
+
+ /* get outputs */
+ E->control.total_iteration_cycles = kE.control.total_iteration_cycles;
+ E->control.total_v_solver_calls = kE.control.total_v_solver_calls;
+ E->monitor.momentum_residual = kE.monitor.momentum_residual;
+
+ free(kE.memory);
+
+ return kE.valid;
}
More information about the CIG-COMMITS
mailing list