[cig-commits] r7761 - in short/3D/PyLith/trunk: .
examples/twocells/twotri3 modulesrc/solver pylith/solver
knepley at geodynamics.org
knepley at geodynamics.org
Sat Jul 28 08:58:33 PDT 2007
Author: knepley
Date: 2007-07-28 08:58:32 -0700 (Sat, 28 Jul 2007)
New Revision: 7761
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/examples/twocells/twotri3/sheardisp.cfg
short/3D/PyLith/trunk/modulesrc/solver/solver.pyxe.src
short/3D/PyLith/trunk/pylith/solver/__init__.py
Log:
Some nonlinear solver work and comments on the TODO list
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-07-28 06:26:01 UTC (rev 7760)
+++ short/3D/PyLith/trunk/TODO 2007-07-28 15:58:32 UTC (rev 7761)
@@ -153,6 +153,8 @@
1. Boundary mesh
+ This works and there is a unit test.
+
2. Construct mesh with higher order cells from mesh with lower order cells.
Many mesh generators do not know how to construct higher order
@@ -165,6 +167,12 @@
Output:
* PETSc Mesh
+ This already works. Our code is written in terms of numBasisFuncs, not
+ vertices. Therefore, if we just use the P2 element from FIAT, our only
+ problem is to correctly allocate the section (I think). To do this, we
+ add the appropriate setFiberDimension() call. FIAT has this information,
+ we just are not using it right now.
+
3. Global refinement of mesh.
Inputs:
@@ -173,6 +181,14 @@
Outputs:
* PETSc Mesh (refined)
+ This already works. You just call
+
+ double maxVolume = 0.01;
+ bool interpolate = false;
+ ALE::Obj<ALE::Mesh> newMesh = ALE::Generator::refineMesh(oldMesh, maxVolume, interpolate);
+
+ There is another version that takes an array of maxVolumes, one fo each cell in oldMesh.
+
======================================================================
QUESTIONS FOR LEIF
======================================================================
Modified: short/3D/PyLith/trunk/examples/twocells/twotri3/sheardisp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/twocells/twotri3/sheardisp.cfg 2007-07-28 06:26:01 UTC (rev 7760)
+++ short/3D/PyLith/trunk/examples/twocells/twotri3/sheardisp.cfg 2007-07-28 15:58:32 UTC (rev 7761)
@@ -20,6 +20,7 @@
# We want an implicit formulation.
formulation = pylith.problems.Implicit
+formulation.solver = pylith.solver.SolverLinear
# This is a container for a boundary condition specified on four sides
# of a rectangular domain.
Modified: short/3D/PyLith/trunk/modulesrc/solver/solver.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/solver/solver.pyxe.src 2007-07-28 06:26:01 UTC (rev 7760)
+++ short/3D/PyLith/trunk/modulesrc/solver/solver.pyxe.src 2007-07-28 15:58:32 UTC (rev 7761)
@@ -291,4 +291,85 @@
return
+# ----------------------------------------------------------------------
+cdef class SolverNonlinear(Solver):
+
+ def __init__(self):
+ """
+ Constructor.
+ """
+ Solver.__init__(self)
+ return
+
+
+ def solve(self, fieldOut, jacobian, fieldIn):
+ """
+ Solve nonlinear system.
+ """
+ # create shim for method 'solve'
+ #embed{ int SolverNonlinear_solve(void* objVptr, void* fieldOutVptr, void* jacobianVptr, void* fieldInVptr, void* scatterVptr, void* vecInVptr, void* vecOutVptr)
+ typedef ALE::Mesh::real_section_type real_section_type;
+ try {
+ assert(0 != objVptr);
+ assert(0 != fieldOutVptr);
+ assert(0 != jacobianVptr);
+ assert(0 != fieldInVptr);
+ assert(0 != scatterVptr);
+ assert(0 != vecInVptr);
+ assert(0 != vecOutVptr);
+
+ KSP* ksp = (KSP*) objVptr;
+ ALE::Obj<real_section_type>* fieldOut =
+ (ALE::Obj<real_section_type>*) fieldOutVptr;
+ Mat* jacobian = (Mat*) jacobianVptr;
+ ALE::Obj<real_section_type>* fieldIn =
+ (ALE::Obj<real_section_type>*) fieldInVptr;
+ VecScatter scatter = (VecScatter) scatterVptr;
+ Vec vecIn = (Vec) vecInVptr;
+ Vec vecOut = (Vec) vecOutVptr;
+
+ PetscErrorCode err = 0;
+ Vec localVec;
+
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF, (*fieldIn)->sizeWithBC(),
+ (*fieldIn)->restrict(), &localVec);CHKERRQ(err);
+ err = VecScatterBegin(scatter, localVec, vecIn, INSERT_VALUES, SCATTER_FORWARD
+ );CHKERRQ(err);
+ err = VecScatterEnd(scatter, localVec, vecIn, INSERT_VALUES, SCATTER_FORWARD
+ ); CHKERRQ(err);
+ err = VecDestroy(localVec); CHKERRQ(err);
+ err = KSPSetOperators(*ksp, *jacobian, *jacobian,
+ DIFFERENT_NONZERO_PATTERN); CHKERRQ(err);
+ err = KSPSolve(*ksp, vecIn, vecOut); CHKERRQ(err);
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF, (*fieldOut)->sizeWithBC(),
+ (*fieldOut)->restrict(), &localVec);CHKERRQ(err);
+ err = VecScatterBegin(scatter, vecOut, localVec, INSERT_VALUES, SCATTER_REVERSE
+ ); CHKERRQ(err);
+ err = VecScatterEnd(scatter, vecOut, localVec, INSERT_VALUES, SCATTER_REVERSE
+ ); CHKERRQ(err);
+ err = VecDestroy(localVec); CHKERRQ(err);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+
+ cdef void* fieldOutVptr
+ cdef void* jacobianVptr
+ cdef void* fieldInVptr
+ fieldOutVptr = PyCObject_AsVoidPtr(fieldOut)
+ jacobianVptr = PyCObject_AsVoidPtr(jacobian)
+ fieldInVptr = PyCObject_AsVoidPtr(fieldIn)
+ SolverNonlinear_solve(self.thisptr, fieldOutVptr, jacobianVptr,
+ fieldInVptr, self.vecScatterVptr,
+ self.vecInVptr, self.vecOutVptr)
+ return
+
+
# End of file
Modified: short/3D/PyLith/trunk/pylith/solver/__init__.py
===================================================================
--- short/3D/PyLith/trunk/pylith/solver/__init__.py 2007-07-28 06:26:01 UTC (rev 7760)
+++ short/3D/PyLith/trunk/pylith/solver/__init__.py 2007-07-28 15:58:32 UTC (rev 7761)
@@ -15,7 +15,8 @@
## @brief Python Pylith solver module initialization
__all__ = ['Solver',
- 'SolverLinear']
+ 'SolverLinear',
+ 'SolverNonlinear']
# End of file
More information about the cig-commits
mailing list