[cig-commits] r14711 - in short/3D/PyLith/branches/pylith-swig: libsrc/problems libsrc/topology modulesrc/problems pylith/problems
knepley at geodynamics.org
knepley at geodynamics.org
Tue Apr 14 19:08:39 PDT 2009
Author: knepley
Date: 2009-04-14 19:08:39 -0700 (Tue, 14 Apr 2009)
New Revision: 14711
Modified:
short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc
short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh
short/3D/PyLith/branches/pylith-swig/libsrc/problems/SolverNonlinear.cc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh
short/3D/PyLith/branches/pylith-swig/modulesrc/problems/Formulation.i
short/3D/PyLith/branches/pylith-swig/pylith/problems/Implicit.py
Log:
Fixing nonlinear solver
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc 2009-04-15 01:31:04 UTC (rev 14710)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc 2009-04-15 02:08:39 UTC (rev 14711)
@@ -87,13 +87,20 @@
_dt = dt;
} // updateSettings
+
// ----------------------------------------------------------------------
// Reform system residual.
void
-pylith::problems::Formulation::reformResidual(void)
+pylith::problems::Formulation::reformResidual(Vec solutionVec, Vec residualVec)
{ // reformResidual
assert(0 != _fields);
+ // Need to pass these Vecs for updating
+
+ // Update section view of field.
+ topology::Field<topology::Mesh>& solution = _fields->get("solution");
+ solution.scatterVectorToSection(solutionVec);
+
// Set residual to zero.
topology::Field<topology::Mesh>& residual = _fields->get("residual");
residual.zero();
@@ -123,17 +130,33 @@
_submeshIntegrators[i]->integrateResidual(residual, _t, _fields);
// Update PETSc view of residual
- residual.scatterSectionToVector();
+ residual.scatterSectionToVector(residualVec);
+
+ // TODO: Move this to SolverLinear
+ VecScale(residualVec, -1.0);
} // reformResidual
// ----------------------------------------------------------------------
// Reform system Jacobian.
void
-pylith::problems::Formulation::reformJacobian()
+pylith::problems::Formulation::reformJacobian(void)
{ // reformJacobian
+ reformJacobian(NULL);
+} // reformJacobian
+void
+pylith::problems::Formulation::reformJacobian(Vec solutionVec)
+{ // reformJacobian
assert(0 != _jacobian);
assert(0 != _fields);
+ // Need to pass these Vecs for updating
+
+ // Update section view of field.
+ if (solutionVec != NULL) {
+ topology::Field<topology::Mesh>& solution = _fields->get("solution");
+ solution.scatterVectorToSection(solutionVec);
+ }
+
// Set residual to zero.
_jacobian->zero();
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh 2009-04-15 01:31:04 UTC (rev 14710)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh 2009-04-15 02:08:39 UTC (rev 14711)
@@ -81,10 +81,11 @@
const double dt);
/// Reform system residual.
- void reformResidual(void);
+ void reformResidual(PetscVec solutionVec, PetscVec residualVec);
/// Reform system Jacobian.
void reformJacobian(void);
+ void reformJacobian(PetscVec solutionVec);
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
@@ -93,6 +94,7 @@
double _dt; ///< Current time step (nondimensional).
topology::Jacobian* _jacobian; ///< Handle to Jacobian of system.
topology::SolutionFields* _fields; ///< Handle to solution fields for system.
+ topology::Field<topology::Mesh>* _solution;
/// Integrators over subdomains of the mesh.
std::vector<IntegratorMesh*> _meshIntegrators;
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/problems/SolverNonlinear.cc 2009-04-15 01:31:04 UTC (rev 14710)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/problems/SolverNonlinear.cc 2009-04-15 02:08:39 UTC (rev 14711)
@@ -93,10 +93,12 @@
const PetscVec residualVec = residual.vector();
const PetscVec solutionVec = solution->vector();
- err = SNESSolve(_snes, residualVec, solutionVec); CHECK_PETSC_ERROR(err);
+ MatView(jacobian.matrix(), PETSC_VIEWER_STDOUT_WORLD);
+ err = SNESSolve(_snes, PETSC_NULL, solutionVec); CHECK_PETSC_ERROR(err);
+ VecView(solutionVec, PETSC_VIEWER_STDOUT_WORLD);
// Update section view of field.
- solution->scatterVectorToSection();
+ solution->scatterVectorToSection(solutionVec);
} // solve
// ----------------------------------------------------------------------
@@ -113,7 +115,8 @@
assert(0 != formulation);
// Reform residual
- formulation->reformResidual();
+ formulation->reformResidual(solutionVec, residualVec);
+ VecView(residualVec, PETSC_VIEWER_STDOUT_WORLD);
return 0;
} // reformResidual
@@ -133,7 +136,7 @@
Formulation* formulation = (Formulation*) context;
assert(0 != formulation);
- formulation->reformJacobian();
+ formulation->reformJacobian(solutionVec);
return 0;
} // reformJacobian
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc 2009-04-15 01:31:04 UTC (rev 14710)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc 2009-04-15 02:08:39 UTC (rev 14711)
@@ -476,18 +476,26 @@
void
pylith::topology::Field<mesh_type>::scatterSectionToVector(void) const
{ // scatterSectionToVector
+ assert(0 != _vector);
+ scatterSectionToVector(_vector);
+} // scatterSectionToVector
+
+template<typename mesh_type>
+void
+pylith::topology::Field<mesh_type>::scatterSectionToVector(Vec vector) const
+{ // scatterSectionToVector
assert(!_section.isNull());
assert(0 != _scatter);
- assert(0 != _vector);
+ assert(0 != vector);
PetscErrorCode err = 0;
PetscVec localVec = 0;
err = VecCreateSeqWithArray(PETSC_COMM_SELF,
_section->sizeWithBC(), _section->restrictSpace(),
&localVec); CHECK_PETSC_ERROR(err);
- err = VecScatterBegin(_scatter, localVec, _vector,
+ err = VecScatterBegin(_scatter, localVec, vector,
INSERT_VALUES, SCATTER_FORWARD); CHECK_PETSC_ERROR(err);
- err = VecScatterEnd(_scatter, localVec, _vector,
+ err = VecScatterEnd(_scatter, localVec, vector,
INSERT_VALUES, SCATTER_FORWARD); CHECK_PETSC_ERROR(err);
err = VecDestroy(localVec); CHECK_PETSC_ERROR(err);
} // scatterSectionToVector
@@ -499,18 +507,25 @@
void
pylith::topology::Field<mesh_type>::scatterVectorToSection(void) const
{ // scatterVectorToSection
+ assert(0 != _vector);
+ scatterVectorToSection(_vector);
+} // scatterVectorToSection
+template<typename mesh_type>
+void
+pylith::topology::Field<mesh_type>::scatterVectorToSection(Vec vector) const
+{ // scatterVectorToSection
assert(!_section.isNull());
assert(0 != _scatter);
- assert(0 != _vector);
+ assert(0 != vector);
PetscErrorCode err = 0;
PetscVec localVec = 0;
err = VecCreateSeqWithArray(PETSC_COMM_SELF,
_section->sizeWithBC(), _section->restrictSpace(),
&localVec); CHECK_PETSC_ERROR(err);
- err = VecScatterBegin(_scatter, _vector, localVec,
+ err = VecScatterBegin(_scatter, vector, localVec,
INSERT_VALUES, SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
- err = VecScatterEnd(_scatter, _vector, localVec,
+ err = VecScatterEnd(_scatter, vector, localVec,
INSERT_VALUES, SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
err = VecDestroy(localVec); CHECK_PETSC_ERROR(err);
} // scatterVectorToSection
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh 2009-04-15 01:31:04 UTC (rev 14710)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh 2009-04-15 02:08:39 UTC (rev 14711)
@@ -229,10 +229,12 @@
/// Scatter section information across processors to update the
/// PETSc vector view of the field.
void scatterSectionToVector(void) const;
+ void scatterSectionToVector(Vec vector) const;
/// Scatter PETSc vector information across processors to update the
/// Sieve section view of the field.
void scatterVectorToSection(void) const;
+ void scatterVectorToSection(Vec vector) const;
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/pylith-swig/modulesrc/problems/Formulation.i
===================================================================
--- short/3D/PyLith/branches/pylith-swig/modulesrc/problems/Formulation.i 2009-04-15 01:31:04 UTC (rev 14710)
+++ short/3D/PyLith/branches/pylith-swig/modulesrc/problems/Formulation.i 2009-04-15 02:08:39 UTC (rev 14711)
@@ -61,7 +61,7 @@
const double dt);
/// Reform system residual.
- void reformResidual(void);
+ void reformResidual(PetscVec solutionVec, PetscVec residualVec);
/// Reform system Jacobian.
void reformJacobian(void);
Modified: short/3D/PyLith/branches/pylith-swig/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/pylith/problems/Implicit.py 2009-04-15 01:31:04 UTC (rev 14710)
+++ short/3D/PyLith/branches/pylith-swig/pylith/problems/Implicit.py 2009-04-15 02:08:39 UTC (rev 14711)
@@ -184,7 +184,7 @@
### NONLINEAR: This moves under SNES control as IntegrateResidual()
### NONLINEAR: Also move updateState() from Integrator.poststep() to this function
- self._reformResidual(t+dt, dt)
+ ### self._reformResidual(t+dt, dt)
self._info.log("Solving equations.")
residual = self.fields.get("residual")
More information about the CIG-COMMITS
mailing list