[cig-commits] r15414 - mc/3D/CitcomS/trunk/lib
leif at geodynamics.org
leif at geodynamics.org
Wed Jul 1 12:09:31 PDT 2009
Author: leif
Date: 2009-07-01 12:09:31 -0700 (Wed, 01 Jul 2009)
New Revision: 15414
Modified:
mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
mc/3D/CitcomS/trunk/lib/cgrad_kernel.cu
mc/3D/CitcomS/trunk/lib/multigrid_kernel.cu
Log:
Fixed a pair of mistakes I made: under cgrad, don't replace
solve_del2_u(); under multigrid, 'level' is not a constant.
Modified: mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2009-07-01 18:51:39 UTC (rev 15413)
+++ mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2009-07-01 19:09:31 UTC (rev 15414)
@@ -45,7 +45,6 @@
Iterative solver also using multigrid ........
=========================================================== */
-#ifndef USE_CUDA
int solve_del2_u(E,d0,F,acc,high_lev)
struct All_variables *E;
double **d0;
@@ -146,7 +145,6 @@
return(valid);
}
-#endif /* !USE_CUDA */
/* =================================
recursive multigrid function ....
Modified: mc/3D/CitcomS/trunk/lib/cgrad_kernel.cu
===================================================================
--- mc/3D/CitcomS/trunk/lib/cgrad_kernel.cu 2009-07-01 18:51:39 UTC (rev 15413)
+++ mc/3D/CitcomS/trunk/lib/cgrad_kernel.cu 2009-07-01 19:09:31 UTC (rev 15414)
@@ -50,14 +50,6 @@
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;
@@ -66,13 +58,6 @@
int n_octoterms;
struct octoterm *ot;
- /* outputs */
-
- struct /*MONITOR*/ {
- double momentum_residual;
- } monitor;
-
- int valid;
};
@@ -367,7 +352,7 @@
/*------------------------------------------------------------------------*/
/* from General_matix_functions.c */
-double conj_grad(
+double do_conj_grad(
struct Some_variables *E,
double *d0,
double *F,
@@ -498,42 +483,6 @@
}
-void do_solve_del2_u(
- struct Some_variables *E,
- double *d0,
- double *F,
- double acc
- )
-{
- int count,cycles;
- int i, neq;
-
- double residual;
-
- neq = E->lmesh.NEQ;
-
- for(i=0;i<neq;i++) {
- d0[i] = 0.0;
- }
-
- residual = sqrt(global_vdot(E,F,F));
-
- count = 0;
-
- cycles = E->control.v_steps_low;
- residual = conj_grad(E,d0,F,acc,&cycles);
- E->valid = (residual < acc)? 1:0;
-
- count++;
-
- E->monitor.momentum_residual = residual;
- E->control.total_iteration_cycles += count;
- E->control.total_v_solver_calls += 1;
-
- return;
-}
-
-
/*------------------------------------------------------------------------*/
static void assert_assumptions(struct All_variables *E, int high_lev) {
@@ -552,17 +501,20 @@
assert(E->parallel.Skip_neq[LEVEL][M] == 0);
}
-extern "C" int solve_del2_u(
+
+extern "C" double conj_grad(
struct All_variables *E,
double **d0,
double **F,
double acc,
- int high_lev
+ int *cycles,
+ int level
)
{
struct Some_variables kE;
+ double residual;
- assert_assumptions(E, high_lev);
+ assert_assumptions(E, level);
/* initialize 'Some_variables' with 'All_variables' */
@@ -584,11 +536,6 @@
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 */
@@ -601,21 +548,12 @@
kE.ot = (struct octoterm *)malloc(kE.n_octoterms * sizeof(struct octoterm));
e_map(&kE);
- /* zero outputs */
- kE.monitor.momentum_residual = 0.0;
- kE.valid = 0;
-
/* solve */
- do_solve_del2_u(&kE, d0[M], F[M], acc);
+ residual = do_conj_grad(&kE, d0[M], F[M], acc, cycles);
- /* 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);
free(kE.mm);
free(kE.ot);
- return kE.valid;
+ return residual;
}
Modified: mc/3D/CitcomS/trunk/lib/multigrid_kernel.cu
===================================================================
--- mc/3D/CitcomS/trunk/lib/multigrid_kernel.cu 2009-07-01 18:51:39 UTC (rev 15413)
+++ mc/3D/CitcomS/trunk/lib/multigrid_kernel.cu 2009-07-01 19:09:31 UTC (rev 15414)
@@ -3,10 +3,11 @@
#include "global_defs.h"
#include "element_definitions.h"
+#include <assert.h>
+#include <stdio.h>
enum {
- LEVEL = 0,
CAPS_PER_PROC = 1,
M = 1, /* cap # */
NSD = 3, /* Spatial extent: 3d */
@@ -224,8 +225,6 @@
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);
}
@@ -246,24 +245,41 @@
/* 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.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.lmesh.NEQ = E->lmesh.NEQ[level];
+ kE.lmesh.NNO = E->lmesh.NNO[level];
- kE.ID = E->ID[LEVEL][M];
+ 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.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.BI = E->BI[level][M];
kE.temp = E->temp[M];
- kE.NODE = E->NODE[LEVEL][M];
+ kE.NODE = E->NODE[level][M];
- /* XXX */
+ /* XXX */
+ do {
+ int i, doff, print;
+ for (i=1;i<=kE.lmesh.NNO;i++) {
+ print = (i < 10 || i > kE.lmesh.NNO - 10);
+ if (print)
+ fprintf(stderr, "%04d:", i);
+ for (doff = 1; doff <= 3; ++doff) {
+ assert(kE.ID[i].doff[doff] == /*NSD*/ 3 * (i - 1) + doff - 1);
+ if (print)
+ fprintf(stderr, " %d", kE.ID[i].doff[doff]);
+ }
+ if (print)
+ fprintf(stderr, "\n");
+ }
+ fprintf(stderr, "\n0 - NEQ %d\n", kE.lmesh.NEQ);
+ } while (0);
+ assert(0);
}
More information about the CIG-COMMITS
mailing list