[cig-commits] r6500 - in short/3D/PyLith/trunk: . modulesrc/solver

[knepley at geodynamics.org]

Author: knepley
Date: 2007-03-31 18:05:41 -0700 (Sat, 31 Mar 2007)
New Revision: 6500

Modified:
   short/3D/PyLith/trunk/TODO
   short/3D/PyLith/trunk/modulesrc/solver/solver.pyxe.src
Log:
Trial code for solving


Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO	2007-04-01 00:51:51 UTC (rev 6499)
+++ short/3D/PyLith/trunk/TODO	2007-04-01 01:05:41 UTC (rev 6500)
@@ -122,7 +122,10 @@
   the appropriate discretization for a problem. For example, I think
   it would allow us to use spectral elements.
 
+6. Trial run with quadratic elements
 
+7. Trial run with fully interpolated mesh
+
 ======================================================================
 QUESTIONS FOR LEIF
 ======================================================================

Modified: short/3D/PyLith/trunk/modulesrc/solver/solver.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/solver/solver.pyxe.src	2007-04-01 00:51:51 UTC (rev 6499)
+++ short/3D/PyLith/trunk/modulesrc/solver/solver.pyxe.src	2007-04-01 01:05:41 UTC (rev 6500)
@@ -60,6 +60,9 @@
     self.name = "pylith_solver_Solver"
 
     # ADD STUFF HERE (CREATE KSP)
+    KSP ksp;
+    ierr = KSPCreate(comm, &ksp);CHKERRQ(ierr);
+    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
     return
 
 
@@ -74,8 +77,10 @@
       assert(0 != meshVptr);
       assert(0 != fieldVptr);
       ALE::Obj<ALE::Field::Mesh>* mesh = (ALE::Obj<ALE::Field::Mesh>*) meshVptr;
-      // STUFF GOES HERE
-      } catch (const std::exception& err) {
+      ALE::Obj<ALE::Field::Mesh::real_section_type>* field = (ALE::Obj<ALE::Field::Mesh::real_section_type>*) fieldVptr;
+      VecScatter scatter;
+      ierr = MeshCreateGlobalScatter(*mesh, *field, &scatter);CHKERRQ(ierr);
+    } catch (const std::exception& err) {
       PyErr_SetString(PyExc_RuntimeError,
                       const_cast<char*>(err.what()));
     } catch (...) {
@@ -123,7 +128,17 @@
       assert(0 != fieldOutVptr);
       assert(0 != jacobianVptr);
       assert(0 != fieldInVptr);
-      // STUFF GOES HERE
+      Vec localVec;
+      ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, fieldIn->size(), fieldIn->restrict(), &localVec);CHKERRQ(ierr);
+      ierr = VecScatterBegin(localVec, vecIn, INSERT_VALUES, SCATTER_FORWARD, scatter);CHKERRQ(ierr);
+      ierr = VecScatterEnd(localVec, vecIn, INSERT_VALUES, SCATTER_FORWARD, scatter);CHKERRQ(ierr);
+      ierr = VecDestroy(localVec);CHKERRQ(ierr);
+      ierr = KSPSetOperators(ksp, jacobian, jacobian, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
+      ierr = KSPSolve(vecIn, vecOut);CHKERRQ(ierr);
+      ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, fieldOut->size(), fieldOut->restrict(), &localVec);CHKERRQ(ierr);
+      ierr = VecScatterBegin(vecOut, localVec, INSERT_VALUES, SCATTER_REVERSE, scatter);CHKERRQ(ierr);
+      ierr = VecScatterEnd(vecOut, localVec, INSERT_VALUES, SCATTER_REVERSE, scatter);CHKERRQ(ierr);
+      ierr = VecDestroy(localVec);CHKERRQ(ierr);
       } catch (const std::exception& err) {
       PyErr_SetString(PyExc_RuntimeError,
                       const_cast<char*>(err.what()));