[cig-commits] r17153 - in mc/3D/CitcomS/trunk: CitcomS/Coupler module/Exchanger

[tan2 at geodynamics.org]

Author: tan2
Date: 2010-08-30 13:41:40 -0700 (Mon, 30 Aug 2010)
New Revision: 17153

Modified:
   mc/3D/CitcomS/trunk/CitcomS/Coupler/Coupler.py
   mc/3D/CitcomS/trunk/module/Exchanger/BaseSVTInlet.cc
   mc/3D/CitcomS/trunk/module/Exchanger/BaseSVTInlet.h
Log:
Added support of multigrid solver in Exchanger. Lifted the solver assertion in Coupler.

Cookbook9 can be converted to multigrid solver once appropriate parameters
(*solver.mesher.levels=2 and *solver.vsolver.Solver=multigrid) are added.
However, the embedded solver will converge incorrectly. Setting
esolver.vsolver.accuracy=6e-3 seems to fix the problem.


Modified: mc/3D/CitcomS/trunk/CitcomS/Coupler/Coupler.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Coupler/Coupler.py	2010-08-30 20:39:06 UTC (rev 17152)
+++ mc/3D/CitcomS/trunk/CitcomS/Coupler/Coupler.py	2010-08-30 20:41:40 UTC (rev 17153)
@@ -55,11 +55,6 @@
         assert solver.inventory.tsolver.inventory.monitor_max_T == False, \
                'Error: solver.tsolver.monitor_max_T must be off!'
 
-        # the exchanger doesn't know how to apply boundary conditions
-        # for the multigrid solver
-        assert solver.inventory.vsolver.inventory.Solver == "cgrad", \
-               'Error: solver.vsolver.Solver must be "cgrad"'
-
         self.communicator = solver.communicator
         self.srcCommList = solver.myPlus
 

Modified: mc/3D/CitcomS/trunk/module/Exchanger/BaseSVTInlet.cc
===================================================================
--- mc/3D/CitcomS/trunk/module/Exchanger/BaseSVTInlet.cc	2010-08-30 20:39:06 UTC (rev 17152)
+++ mc/3D/CitcomS/trunk/module/Exchanger/BaseSVTInlet.cc	2010-08-30 20:41:40 UTC (rev 17153)
@@ -39,7 +39,7 @@
 #include "element_definitions.h"
 
     void check_bc_consistency(const All_variables *E);
-    void construct_id(const All_variables *E);
+    void get_bcs_id_for_residual(const All_variables *E, int level, int m);
 }
 
 using Exchanger::Array2D;
@@ -103,18 +103,24 @@
 
     for(int i=0; i<boundary.size(); i++) {
         int n = boundary.nodeID(i);
-	for(int d=0; d<DIM; ++d)
-	    if(boundary.normal(d,i)) {
-		E->node[m][n] = E->node[m][n] | vbcFlag[d];
-		E->node[m][n] = E->node[m][n] & (~sbcFlag[d]);
-	    } else {
-		E->node[m][n] = E->node[m][n] | sbcFlag[d];
-		E->node[m][n] = E->node[m][n] & (~vbcFlag[d]);
-	    }
+        for(int lev=E->mesh.gridmax; lev>=E->mesh.gridmin; lev--) {
+            int nn = nodelevel(n, lev);
+            if(nn != 0) {
+                for(int d=0; d<DIM; ++d)
+                    if(boundary.normal(d,i)) {
+                        E->NODE[lev][m][nn] = E->NODE[lev][m][nn] | vbcFlag[d];
+                        E->NODE[lev][m][nn] = E->NODE[lev][m][nn] & (~sbcFlag[d]);
+                    } else {
+                        E->NODE[lev][m][nn] = E->NODE[lev][m][nn] | sbcFlag[d];
+                        E->NODE[lev][m][nn] = E->NODE[lev][m][nn] & (~vbcFlag[d]);
+                    }
+            }
+        }
     }
 
-    // reconstruct ID array to reflect changes in VBC
-    construct_id(E);
+    // reconstruct zero_resid array to reflect changes in VBC
+    for(int lev=E->mesh.gridmax; lev>=E->mesh.gridmin; lev--)
+        get_bcs_id_for_residual(E,lev,m);
 }
 
 
@@ -213,6 +219,34 @@
 }
 
 
+int BaseSVTInlet::nodelevel(int node, int level)
+{
+    // Given the "node" number of highest level, return the node number
+    // at "level" below. Return 0 if this node is not below to lower level.
+    const int nox = E->lmesh.nox;
+    const int noy = E->lmesh.noy;
+    const int noz = E->lmesh.noz;
+    const int nx = E->lmesh.NOX[level];
+    const int ny = E->lmesh.NOY[level];
+    const int nz = E->lmesh.NOZ[level];
+    const int f = pow(2, (E->mesh.gridmax - level));
+
+    int i, j, k, ii, jj, kk;
+
+    // indices at highest level, starting from 0
+    k = (node - 1) % noz;
+    j = (node - 1) / (nox * noz);
+    i = ((node - 1) - j * nox * noz) / noz;
+
+    if((i % f) != 0 || (j % f) != 0 || (k % f) != 0) return 0;
+
+    ii = i / f;
+    jj = j / f;
+    kk = k / f;
+
+    return jj*nx*nz + ii*nz + kk + 1;
+}
+
 // version
 // $Id: BaseSVTInlet.cc 7643 2007-07-11 20:17:32Z tan2 $
 

Modified: mc/3D/CitcomS/trunk/module/Exchanger/BaseSVTInlet.h
===================================================================
--- mc/3D/CitcomS/trunk/module/Exchanger/BaseSVTInlet.h	2010-08-30 20:39:06 UTC (rev 17152)
+++ mc/3D/CitcomS/trunk/module/Exchanger/BaseSVTInlet.h	2010-08-30 20:41:40 UTC (rev 17153)
@@ -71,6 +71,8 @@
 
     double side_tractions(const Exchanger::Array2D<double,Exchanger::STRESS_DIM>& stress,
 			  int node, int normal_dir,  int dim) const;
+
+    int nodelevel(int node, int level);
 };