[cig-commits] r7329 - short/3D/PyLith/trunk/pylith/problems

[brad at geodynamics.org]

Author: brad
Date: 2007-06-20 16:19:15 -0700 (Wed, 20 Jun 2007)
New Revision: 7329

Modified:
   short/3D/PyLith/trunk/pylith/problems/Explicit.py
   short/3D/PyLith/trunk/pylith/problems/Implicit.py
   short/3D/PyLith/trunk/pylith/problems/TimeDependent.py
Log:
Fixed time stepping start time. For implicit time stepping first time step is -dt to 0; for explicit time stepping first time step is 0 to dt. First solution written in each case should be at time 0.

Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py	2007-06-20 22:49:32 UTC (rev 7328)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py	2007-06-20 23:19:15 UTC (rev 7329)
@@ -86,6 +86,14 @@
     return
 
 
+  def startTime(self):
+    """
+    Get time at which time stepping should start.
+    """
+    from pyre.units.time import second
+    return 0.0*second
+
+
   def stableTimeStep(self):
     """
     Get stable time step for advancing forward in time.

Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py	2007-06-20 22:49:32 UTC (rev 7328)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py	2007-06-20 23:19:15 UTC (rev 7329)
@@ -104,6 +104,13 @@
     return
 
 
+  def startTime(self, dt):
+    """
+    Get time at which time stepping should start.
+    """
+    return -dt
+
+
   def stableTimeStep(self):
     """
     Get stable time step for advancing forward in time.
@@ -155,10 +162,9 @@
       integrator.timeStep(dt)
       integrator.integrateResidual(residual, t, self.fields)
 
-    bindings.sectionView(residual, "residual")
+    import pylith.utils.petsc as petsc
     self._info.log("Solving equations.")
     self.solver.solve(dispIncr, self.jacobian, residual)
-    bindings.sectionView(self.fields.getReal("solution"), "solution")
     return
 
 
@@ -175,7 +181,6 @@
     solution = self.fields.getReal("solution")
     disp = self.fields.getReal("dispTBctpdt")
     bindings.addRealSections(disp, disp, solution)
-    bindings.sectionView(solution, "solution")
 
     self._info.log("Updating integrators states.")
     for integrator in self.integrators:
@@ -205,54 +210,6 @@
     return
 
 
-  def _solveElastic(self, mesh, materials, t, dt):
-    """
-    Solve for elastic solution.
-    """
-    self._info.log("Computing elastic solution.")
-
-    self._info.log("Setting constraints.")
-    solnField = self.fields.getReal("dispTBctpdt")
-    import pylith.topology.topology as bindings
-    bindings.zeroRealSection(solnField)
-    for constraint in self.constraints:
-      constraint.setField(t, solnField)
-
-    self._info.log("Integrating Jacobian and residual of operator.")
-    import pylith.utils.petsc as petsc
-    petsc.mat_setzero(self.jacobian)
-    residual = self.fields.getReal("residual")
-    petsc.mat_setzero(self.jacobian)
-    bindings.zeroRealSection(residual)
-    for integrator in self.integrators:
-      integrator.timeStep(dt)
-      integrator.integrateJacobian(self.jacobian, t, self.fields)
-      integrator.integrateResidual(residual, t, self.fields)
-    import pylith.utils.petsc as petsc
-    petsc.mat_assemble(self.jacobian)
-
-    import pylith.topology.topology as bindings
-    petsc.mat_view(self.jacobian)
-    bindings.sectionView(residual, "residual")
-
-    self._info.log("Solving equations.")
-    self.solver.solve(solution, self.jacobian, residual)
-
-    bindings.addRealSections(disp, disp, solution)
-    bindings.sectionView(solution, "solution")
-    bindings.sectionView(self.fields.getReal("dispIncr"), "dispIncr")
-
-    self._info.log("Updating integrators states.")
-    for integrator in self.integrators:
-      integrator.updateState(t, disp)
-
-    self._info.log("Outputting elastic solution.")
-    for output in self.output.bin:
-      output.writeField(t, self._istep, disp, self.solnField['label'])
-    self._istep += 1      
-    return
-  
-
 # FACTORIES ////////////////////////////////////////////////////////////
 
 def pde_formulation():

Modified: short/3D/PyLith/trunk/pylith/problems/TimeDependent.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/TimeDependent.py	2007-06-20 22:49:32 UTC (rev 7328)
+++ short/3D/PyLith/trunk/pylith/problems/TimeDependent.py	2007-06-20 23:19:15 UTC (rev 7329)
@@ -109,10 +109,11 @@
     self._info.log("Solving problem.")
     self.checkpointTimer.toplevel = app # Set handle for saving state
     
-    from pyre.units.time import second
-    t = 0.0*second
-
-    while t.value <= self.totalTime.value:
+    dt = self.formulation.stableTimeStep()
+    if dt.value == 0.0: # If formulation returns 0.0, use default time step
+      dt = self.dt
+    t = self.formulation.startTime(self.dt)
+    while t.value < self.totalTime.value:
       self._info.log("Main time loop, t=%s" % t)
       
       # Checkpoint if necessary