[cig-commits] r19813 - in short/3D/PyLith/branches/pylith-scecdynrup: applications applications/utilities examples/3d/hex8 libsrc/pylith/problems libsrc/pylith/utils
brad at geodynamics.org
brad at geodynamics.org
Mon Mar 19 11:15:38 PDT 2012
Author: brad
Date: 2012-03-19 11:15:38 -0700 (Mon, 19 Mar 2012)
New Revision: 19813
Added:
short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylith_genxdmf.in
short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylithinfo.in
Removed:
short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylith_genxdmf
short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylithinfo
Modified:
short/3D/PyLith/branches/pylith-scecdynrup/applications/pylith.in
short/3D/PyLith/branches/pylith-scecdynrup/applications/pylith_prepmesh.in
short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.cfg
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.hh
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/utils/petscfwd.h
Log:
Merge from trunk.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/applications/pylith.in
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/applications/pylith.in 2012-03-19 18:14:13 UTC (rev 19812)
+++ short/3D/PyLith/branches/pylith-scecdynrup/applications/pylith.in 2012-03-19 18:15:38 UTC (rev 19813)
@@ -10,7 +10,7 @@
# This code was developed as part of the Computational Infrastructure
# for Geodynamics (http://geodynamics.org).
#
-# Copyright (c) 2010-2011 University of California, Davis
+# Copyright (c) 2010-2012 University of California, Davis
#
# See COPYING for license information.
#
Modified: short/3D/PyLith/branches/pylith-scecdynrup/applications/pylith_prepmesh.in
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/applications/pylith_prepmesh.in 2012-03-19 18:14:13 UTC (rev 19812)
+++ short/3D/PyLith/branches/pylith-scecdynrup/applications/pylith_prepmesh.in 2012-03-19 18:15:38 UTC (rev 19813)
@@ -10,7 +10,7 @@
# This code was developed as part of the Computational Infrastructure
# for Geodynamics (http://geodynamics.org).
#
-# Copyright (c) 2010-2011 University of California, Davis
+# Copyright (c) 2010-2012 University of California, Davis
#
# See COPYING for license information.
#
Modified: short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/Makefile.am 2012-03-19 18:14:13 UTC (rev 19812)
+++ short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/Makefile.am 2012-03-19 18:15:38 UTC (rev 19813)
@@ -16,10 +16,50 @@
# ----------------------------------------------------------------------
#
-dist_bin_SCRIPTS = \
+bin_SCRIPTS = \
pylithinfo \
pylith_genxdmf \
powerlaw_gendb.py
+INTERPRETER = $(PYTHON)
+noinstINTERPRETER = $(PYTHON)
+
+do_build = sed -e s%[@]INTERPRETER[@]%$(noinstINTERPRETER)%g -e s%[@]PYTHONPATH[@]%$(noinstPYTHONPATH)%g
+do_install = sed -e s%[@]INTERPRETER[@]%$(INTERPRETER)%g -e s%[@]PYTHONPATH[@]%$(PYTHONPATH)%g
+
+pylithinfo: $(srcdir)/pylithinfo.in Makefile
+ $(do_build) < $(srcdir)/pylithinfo.in > $@ || (rm -f $@ && exit 1)
+ chmod +x $@
+
+pylith_genxdmf: $(srcdir)/pylith_genxdmf.in Makefile
+ $(do_build) < $(srcdir)/pylith_genxdmf.in > $@ || (rm -f $@ && exit 1)
+ chmod +x $@
+
+install-binSCRIPTS: $(bin_SCRIPTS)
+ @$(NORMAL_INSTALL)
+ test -z "$(bindir)" || $(mkdir_p) "$(DESTDIR)$(bindir)"
+ @list='$(bin_SCRIPTS)'; for p in $$list; do \
+ if test -f "$$p.in"; then d=; else d="$(srcdir)/"; fi; \
+ if test -f $$d$$p.in; then \
+ f=`echo "$$p" | sed 's|^.*/||;$(transform)'`; \
+ echo " $(do_install) '$$d$$p.in' > '$(DESTDIR)$(bindir)/$$f'"; \
+ $(do_install) "$$d$$p.in" > "$(DESTDIR)$(bindir)/$$f"; \
+ echo " chmod +x '$(DESTDIR)$(bindir)/$$f'"; \
+ chmod +x "$(DESTDIR)$(bindir)/$$f"; \
+ else :; fi; \
+ done
+
+EXTRA_DIST = \
+ pylithinfo.in \
+ pylith_genxdmf.in
+
+CLEANFILES = $(bin_SCRIPTS)
+
+
# End of file
+
+
+
+
+# End of file
Deleted: short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylith_genxdmf
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylith_genxdmf 2012-03-19 18:14:13 UTC (rev 19812)
+++ short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylith_genxdmf 2012-03-19 18:15:38 UTC (rev 19813)
@@ -1,74 +0,0 @@
-#!/usr/bin/env nemesis
-#
-# ======================================================================
-#
-# Brad T. Aagaard, U.S. Geological Survey
-# Charles A. Williams, GNS Science
-# Matthew G. Knepley, University of Chicago
-#
-# This code was developed as part of the Computational Infrastructure
-# for Geodynamics (http://geodynamics.org).
-#
-# Copyright (c) 2010-2011 University of California, Davis
-#
-# See COPYING for license information.
-#
-# ======================================================================
-#
-
-# This script creates a metadata file (.xmf) from an HDF5 file written
-# by PyLith permitting the data to be viewed in VTK visualization
-# software, such as ParaView and Visit.
-#
-# Usage: pylith_genxdmf --file=FILE
-
-__requires__ = "PyLith"
-
-
-# ======================================================================
-class GenXdmfApp(object):
- """
- Application for generating an Xdmf file for an HDF5 file.
- """
-
- # PUBLIC METHODS /////////////////////////////////////////////////////
-
- def __init__(self, name="pylith_genxdmf"):
- """
- Constructor.
- """
- self.filename = "output.h5"
- return
-
-
- def main(self, *args, **kwds):
- """
- Main entry point for application.
- """
- from pylith.meshio.Xdmf import Xdmf
- writer = Xdmf()
- filenameHDF5 = self.filename
- filenameXdmf = filenameHDF5.replace(".h5", ".xmf")
- writer.write(filenameXdmf, filenameHDF5)
- return
-
-
-# ----------------------------------------------------------------------
-if __name__ == "__main__":
-
- usage = "%prog --file=FILE"
- from optparse import OptionParser
- parser = OptionParser(usage=usage)
- parser.add_option("-f", "--file",
- dest="filename",
- type="string", metavar="FILE",
- help="Create Xdmf file for HDF5 file FILE. [output.h5]",
- default="output.h5")
- (options, args) = parser.parse_args()
-
- app = GenXdmfApp()
- app.filename = options.filename
- app.main()
-
-
-# End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylith_genxdmf.in (from rev 19812, short/3D/PyLith/trunk/applications/utilities/pylith_genxdmf.in)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylith_genxdmf.in (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylith_genxdmf.in 2012-03-19 18:15:38 UTC (rev 19813)
@@ -0,0 +1,75 @@
+#!@INTERPRETER@
+# -*- Python -*-
+#
+# ======================================================================
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ======================================================================
+#
+
+# This script creates a metadata file (.xmf) from an HDF5 file written
+# by PyLith permitting the data to be viewed in VTK visualization
+# software, such as ParaView and Visit.
+#
+# Usage: pylith_genxdmf --file=FILE
+
+__requires__ = "PyLith"
+
+
+# ======================================================================
+class GenXdmfApp(object):
+ """
+ Application for generating an Xdmf file for an HDF5 file.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="pylith_genxdmf"):
+ """
+ Constructor.
+ """
+ self.filename = "output.h5"
+ return
+
+
+ def main(self, *args, **kwds):
+ """
+ Main entry point for application.
+ """
+ from pylith.meshio.Xdmf import Xdmf
+ writer = Xdmf()
+ filenameHDF5 = self.filename
+ filenameXdmf = filenameHDF5.replace(".h5", ".xmf")
+ writer.write(filenameXdmf, filenameHDF5)
+ return
+
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+
+ usage = "%prog --file=FILE"
+ from optparse import OptionParser
+ parser = OptionParser(usage=usage)
+ parser.add_option("-f", "--file",
+ dest="filename",
+ type="string", metavar="FILE",
+ help="Create Xdmf file for HDF5 file FILE. [output.h5]",
+ default="output.h5")
+ (options, args) = parser.parse_args()
+
+ app = GenXdmfApp()
+ app.filename = options.filename
+ app.main()
+
+
+# End of file
Deleted: short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylithinfo
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylithinfo 2012-03-19 18:14:13 UTC (rev 19812)
+++ short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylithinfo 2012-03-19 18:15:38 UTC (rev 19813)
@@ -1,193 +0,0 @@
-#!/usr/bin/env nemesis
-#
-# ======================================================================
-#
-# Brad T. Aagaard, U.S. Geological Survey
-# Charles A. Williams, GNS Science
-# Matthew G. Knepley, University of Chicago
-#
-# This code was developed as part of the Computational Infrastructure
-# for Geodynamics (http://geodynamics.org).
-#
-# Copyright (c) 2010-2011 University of California, Davis
-#
-# See COPYING for license information.
-#
-# ======================================================================
-#
-
-# This script dumps all PyLith parameters (defaults plus those
-# specified by the user to a text file. The default name of the output
-# file is 'pylith_parameters.txt'. Verbose output includes a
-# description of the parameter along with where it's current value was
-# set.
-#
-# Usage: pylithinfo.py [--verbose] [--fileout=FILE] [PyLith args]
-
-__requires__ = "PyLith"
-
-
-# ======================================================================
-class ParametersApp(object):
- """
- Application for printing current PyLith parameters to a text file.
- """
-
- # PUBLIC METHODS /////////////////////////////////////////////////////
-
- def __init__(self, name="pylithinfo"):
- """
- Constructor.
- """
- self.verbose = False
- self.filename = "pylith_parameters.txt"
- self.pylith_args = ""
- self._tab = " "
- self._printPropertyFn = self._printPropertyBasic
- self._printFacilityFn = self._printFacilityBasic
- return
-
-
- def main(self, *args, **kwds):
- """
- Main entry point for application.
- """
- from pylith.apps.PyLithApp import InfoApp
- targetapp = InfoApp(self.pylith_args)
- targetapp.run(*args, **kwds)
-
- if self.verbose:
- self._printPropertyFn = self._printPropertyVerbose
- self._printFacilityFn = self._printFacilityVerbose
- else:
- self._printPropertyFn = self._printPropertyBasic
- self._printFacilityFn = self._printFacilityBasic
-
- depth = 0
- fout = open(self.filename, "w")
- fout.write("Application: %s %s\n" % (targetapp.name, targetapp))
- self._printParams(fout, targetapp, depth+1)
- fout.close()
- return
-
-
- # PRIVATE METHODS ////////////////////////////////////////////////////
-
- def _printParams(self, fout, obj, depth):
- """
- Print objects parameters to fout.
- """
- propertyNames = obj.inventory.propertyNames()
- facilityNames = obj.inventory.facilityNames()
-
- propertiesOmit = ["help",
- "help-components",
- "help-persistence",
- "help-properties",
- "typos",
- ]
- facilitiesOmit = ["weaver",
- ]
-
- propertyNames.sort()
- facilityNames.sort()
-
- for name in propertyNames:
- if name in facilityNames or name in propertiesOmit:
- continue
- trait = obj.inventory.getTrait(name)
- descriptor = obj.inventory.getTraitDescriptor(name)
- self._printPropertyFn(fout, name, trait, descriptor, depth)
- for name in facilityNames:
- if name in facilitiesOmit:
- continue
- trait = obj.inventory.getTrait(name)
- descriptor = obj.inventory.getTraitDescriptor(name)
- self._printFacilityFn(fout, name, trait, descriptor, depth)
- return
-
-
- def _printPropertyBasic(self, fout, name, trait, descriptor, depth):
- """
- Print property name, type, and value.
- """
- indent = self._tab*depth
- fout.write("%s%s (%s) = %s\n" % \
- (indent, name, trait.type, descriptor.value))
- return
-
-
- def _printFacilityBasic(self, fout, name, trait, descriptor, depth):
- """
- Print facility name, type, and value.
- """
- indent = self._tab*depth
- fout.write("%s%s = %s (%s)\n" % \
- (indent, name, descriptor.value.name, descriptor.value))
- self._printParams(fout, descriptor.value, depth+1)
- return
-
-
- def _printPropertyVerbose(self, fout, name, trait, descriptor, depth):
- """
- Print property, name, type, value, description, and location set.
- """
- indent = self._tab*depth
- fout.write("\n%s%s (%s) = %s\n" % \
- (indent, name, trait.type, descriptor.value))
-
- indent += self._tab
- try:
- description = trait.meta['tip']
- except KeyError:
- description = "No description available."
- fout.write("%sDescription: %s\n" % (indent, description))
- fout.write("%sSet from: %s\n" % (indent, descriptor.locator))
- return
-
-
- def _printFacilityVerbose(self, fout, name, trait, descriptor, depth):
- """
- Print facility name, type, value, description, and location set.
- """
- indent = self._tab*depth
- fout.write("\n%s%s = %s (%s)\n" % \
- (indent, name, descriptor.value.name, descriptor.value))
-
- indent += self._tab
- try:
- description = trait.meta['tip']
- except KeyError:
- description = "No description available."
- fout.write("%sDescription: %s\n" % (indent, description))
- fout.write("%sSet from: %s\n" % (indent, descriptor.locator))
- fout.write("%sConfigurable as: %s\n" % \
- (indent, ", ".join(descriptor.value.aliases)))
-
- self._printParams(fout, descriptor.value, depth+1)
- return
-
-
-# ----------------------------------------------------------------------
-if __name__ == "__main__":
-
- usage = "%prog [options] [PyLith args]"
- from optparse import OptionParser
- parser = OptionParser(usage=usage)
- parser.add_option("-v", "--verbose", dest="verbose",
- action="store_true", default=False,
- help="Print verbose output.")
- parser.add_option("-o", "--fileout", dest="filename",
- type="string", metavar="FILE",
- help="Write to FILE. [pylith_parameters.txt]",
- default="pylith_parameters.txt")
- (options, args) = parser.parse_args()
-
- app = ParametersApp()
- app.verbose = options.verbose
- app.filename = options.filename
- app.pylith_args = args
- app.main()
-
-
-# End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylithinfo.in (from rev 19812, short/3D/PyLith/trunk/applications/utilities/pylithinfo.in)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylithinfo.in (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/applications/utilities/pylithinfo.in 2012-03-19 18:15:38 UTC (rev 19813)
@@ -0,0 +1,194 @@
+#!@INTERPRETER@
+# -*- Python -*-
+#
+# ======================================================================
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ======================================================================
+#
+
+# This script dumps all PyLith parameters (defaults plus those
+# specified by the user to a text file. The default name of the output
+# file is 'pylith_parameters.txt'. Verbose output includes a
+# description of the parameter along with where it's current value was
+# set.
+#
+# Usage: pylithinfo.py [--verbose] [--fileout=FILE] [PyLith args]
+
+__requires__ = "PyLith"
+
+
+# ======================================================================
+class ParametersApp(object):
+ """
+ Application for printing current PyLith parameters to a text file.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="pylithinfo"):
+ """
+ Constructor.
+ """
+ self.verbose = False
+ self.filename = "pylith_parameters.txt"
+ self.pylith_args = ""
+ self._tab = " "
+ self._printPropertyFn = self._printPropertyBasic
+ self._printFacilityFn = self._printFacilityBasic
+ return
+
+
+ def main(self, *args, **kwds):
+ """
+ Main entry point for application.
+ """
+ from pylith.apps.PyLithApp import InfoApp
+ targetapp = InfoApp(self.pylith_args)
+ targetapp.run(*args, **kwds)
+
+ if self.verbose:
+ self._printPropertyFn = self._printPropertyVerbose
+ self._printFacilityFn = self._printFacilityVerbose
+ else:
+ self._printPropertyFn = self._printPropertyBasic
+ self._printFacilityFn = self._printFacilityBasic
+
+ depth = 0
+ fout = open(self.filename, "w")
+ fout.write("Application: %s %s\n" % (targetapp.name, targetapp))
+ self._printParams(fout, targetapp, depth+1)
+ fout.close()
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _printParams(self, fout, obj, depth):
+ """
+ Print objects parameters to fout.
+ """
+ propertyNames = obj.inventory.propertyNames()
+ facilityNames = obj.inventory.facilityNames()
+
+ propertiesOmit = ["help",
+ "help-components",
+ "help-persistence",
+ "help-properties",
+ "typos",
+ ]
+ facilitiesOmit = ["weaver",
+ ]
+
+ propertyNames.sort()
+ facilityNames.sort()
+
+ for name in propertyNames:
+ if name in facilityNames or name in propertiesOmit:
+ continue
+ trait = obj.inventory.getTrait(name)
+ descriptor = obj.inventory.getTraitDescriptor(name)
+ self._printPropertyFn(fout, name, trait, descriptor, depth)
+ for name in facilityNames:
+ if name in facilitiesOmit:
+ continue
+ trait = obj.inventory.getTrait(name)
+ descriptor = obj.inventory.getTraitDescriptor(name)
+ self._printFacilityFn(fout, name, trait, descriptor, depth)
+ return
+
+
+ def _printPropertyBasic(self, fout, name, trait, descriptor, depth):
+ """
+ Print property name, type, and value.
+ """
+ indent = self._tab*depth
+ fout.write("%s%s (%s) = %s\n" % \
+ (indent, name, trait.type, descriptor.value))
+ return
+
+
+ def _printFacilityBasic(self, fout, name, trait, descriptor, depth):
+ """
+ Print facility name, type, and value.
+ """
+ indent = self._tab*depth
+ fout.write("%s%s = %s (%s)\n" % \
+ (indent, name, descriptor.value.name, descriptor.value))
+ self._printParams(fout, descriptor.value, depth+1)
+ return
+
+
+ def _printPropertyVerbose(self, fout, name, trait, descriptor, depth):
+ """
+ Print property, name, type, value, description, and location set.
+ """
+ indent = self._tab*depth
+ fout.write("\n%s%s (%s) = %s\n" % \
+ (indent, name, trait.type, descriptor.value))
+
+ indent += self._tab
+ try:
+ description = trait.meta['tip']
+ except KeyError:
+ description = "No description available."
+ fout.write("%sDescription: %s\n" % (indent, description))
+ fout.write("%sSet from: %s\n" % (indent, descriptor.locator))
+ return
+
+
+ def _printFacilityVerbose(self, fout, name, trait, descriptor, depth):
+ """
+ Print facility name, type, value, description, and location set.
+ """
+ indent = self._tab*depth
+ fout.write("\n%s%s = %s (%s)\n" % \
+ (indent, name, descriptor.value.name, descriptor.value))
+
+ indent += self._tab
+ try:
+ description = trait.meta['tip']
+ except KeyError:
+ description = "No description available."
+ fout.write("%sDescription: %s\n" % (indent, description))
+ fout.write("%sSet from: %s\n" % (indent, descriptor.locator))
+ fout.write("%sConfigurable as: %s\n" % \
+ (indent, ", ".join(descriptor.value.aliases)))
+
+ self._printParams(fout, descriptor.value, depth+1)
+ return
+
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+
+ usage = "%prog [options] [PyLith args]"
+ from optparse import OptionParser
+ parser = OptionParser(usage=usage)
+ parser.add_option("-v", "--verbose", dest="verbose",
+ action="store_true", default=False,
+ help="Print verbose output.")
+ parser.add_option("-o", "--fileout", dest="filename",
+ type="string", metavar="FILE",
+ help="Write to FILE. [pylith_parameters.txt]",
+ default="pylith_parameters.txt")
+ (options, args) = parser.parse_args()
+
+ app = ParametersApp()
+ app.verbose = options.verbose
+ app.filename = options.filename
+ app.pylith_args = args
+ app.main()
+
+
+# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.cfg 2012-03-19 18:14:13 UTC (rev 19812)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.cfg 2012-03-19 18:15:38 UTC (rev 19813)
@@ -103,7 +103,7 @@
snes_atol = 1.0e-9
snes_max_it = 100
snes_monitor = true
-snes_ls_monitor = true
+linesearch_monitor = true
#snes_view = true
snes_converged_reason = true
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc 2012-03-19 18:14:13 UTC (rev 19812)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc 2012-03-19 18:15:38 UTC (rev 19813)
@@ -35,8 +35,23 @@
#define isnan std::isnan // TEMPORARY
#define isinf std::isinf // TEMPORARY
+// Customized line search based on PETSc code in
+// src/snes/linesearch/bt/linesearchbt.c.
#include <private/snesimpl.h>
+#include <private/linesearchimpl.h>
+typedef enum {
+ SNES_LINESEARCH_BT_QUADRATIC,
+ SNES_LINESEARCH_BT_CUBIC,
+} PetscSNESLineSearchBTOrder;
+
+struct PetscSNESLineSearch_BT {
+ PetscReal alpha; /* sufficient decrease parameter */
+ PetscSNESLineSearchBTOrder order;
+};
+
+#define PYLITH_CUSTOM_LINESEARCH 1
+
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
@@ -95,9 +110,13 @@
CHECK_PETSC_ERROR(err);
err = SNESSetFromOptions(_snes); CHECK_PETSC_ERROR(err);
- err = SNESLineSearchSet(_snes, lineSearch, (void*) formulation); CHECK_PETSC_ERROR(err);
err = SNESSetComputeInitialGuess(_snes, initialGuess, (void*) formulation); CHECK_PETSC_ERROR(err);
+ PetscSNESLineSearch ls;
+ err = SNESGetSNESLineSearch(_snes, &ls); CHECK_PETSC_ERROR(err);
+ err = SNESLineSearchSetType(ls, SNES_LINESEARCH_SHELL); CHECK_PETSC_ERROR(err);
+ err = SNESLineSearchShellSetUserFunc(ls, lineSearch, (void*) formulation); CHECK_PETSC_ERROR(err);
+
if (formulation->splitFields()) {
PetscKSP ksp = 0;
PetscPC pc = 0;
@@ -185,246 +204,297 @@
#undef __FUNCT__
#define __FUNCT__ "lineSearch"
PetscErrorCode
-pylith::problems::SolverNonlinear::lineSearch(PetscSNES snes,
- void *lsctx,
- PetscVec x,
- PetscVec f,
- PetscVec g,
- PetscVec y,
- PetscVec w,
- PetscReal fnorm,
- PetscReal xnorm,
- PetscReal *ynorm,
- PetscReal *gnorm,
- PetscBool *flag)
+pylith::problems::SolverNonlinear::lineSearch(PetscSNESLineSearch linesearch,
+ void* lsctx)
{ // lineSearch
// Note that for line search purposes we work with with the related
// minimization problem:
// min z(x): R^n -> R,
// where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2.
+ //
+ // Customized version of src/snes/linesearch/impls/bt/linesearchbt.c.
- PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength;
- PetscReal minlambda,lambda,lambdatemp;
+ PetscBool changed_y =PETSC_FALSE, changed_w =PETSC_FALSE;
+ PetscErrorCode ierr;
+ Vec X, F, Y, W, G;
+ SNES snes;
+ PetscReal fnorm, xnorm, ynorm, gnorm, gnormprev;
+ PetscReal lambda, lambdatemp, lambdaprev, minlambda, maxstep, rellength, initslope, alpha;
+ PetscReal t1, t2, a, b, d, steptol;
#if defined(PETSC_USE_COMPLEX)
PetscScalar cinitslope;
#endif
- PetscErrorCode ierr;
- PetscInt count;
- PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE;
- MPI_Comm comm;
+ PetscBool domainerror;
+ PetscViewer monitor;
+ PetscInt max_its, count;
+ PetscSNESLineSearch_BT *bt;
+ Mat jac;
+
PetscFunctionBegin;
- ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr);
- ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
- *flag = PETSC_TRUE;
- ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr);
- if (*ynorm == 0.0) {
- if (snes->ls_monitor) {
- ierr = PetscViewerASCIIAddTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
- ierr = PetscViewerASCIIPrintf(snes->ls_monitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr);
- ierr = PetscViewerASCIISubtractTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
+ ierr = SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);CHKERRQ(ierr);
+ ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr);
+ ierr = SNESLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr);
+ ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
+ ierr = SNESLineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr);
+ ierr = SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, PETSC_NULL, PETSC_NULL, PETSC_NULL, &max_its);
+ bt = (PetscSNESLineSearch_BT *)linesearch->data;
+
+ alpha = bt->alpha;
+
+ ierr = SNESGetJacobian(snes, &jac, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr);
+ if (!jac) {
+ SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
+ }
+ /* precheck */
+ ierr = SNESLineSearchPreCheck(linesearch, &changed_y);CHKERRQ(ierr);
+ ierr = SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr);
+
+ ierr = VecNorm(Y, NORM_2, &ynorm);CHKERRQ(ierr);
+ if (ynorm == 0.0) {
+ if (monitor) {
+ ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ ierr = PetscViewerASCIIPrintf(monitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr);
+ ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
}
- *gnorm = fnorm;
- ierr = VecCopy(x,w);CHKERRQ(ierr);
- ierr = VecCopy(f,g);CHKERRQ(ierr);
- *flag = PETSC_FALSE;
- goto theend1;
+ ierr = VecCopy(X,W);CHKERRQ(ierr);
+ ierr = VecCopy(F,G);CHKERRQ(ierr);
+ ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
}
- if (*ynorm > snes->maxstep) { /* Step too big, so scale back */
- if (snes->ls_monitor) {
- ierr = PetscViewerASCIIAddTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
- ierr = PetscViewerASCIIPrintf(snes->ls_monitor," Line search: Scaling step by %14.12e old ynorm %14.12e\n",(double)(snes->maxstep/(*ynorm)),(double)(*ynorm));CHKERRQ(ierr);
- ierr = PetscViewerASCIISubtractTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
+ if (ynorm > maxstep) { /* Step too big, so scale back */
+ if (monitor) {
+ ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ ierr = PetscViewerASCIIPrintf(monitor," Line search: Scaling step by %14.12e old ynorm %14.12e\n", (maxstep/ynorm),ynorm);CHKERRQ(ierr);
+ ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
}
- ierr = VecScale(y,snes->maxstep/(*ynorm));CHKERRQ(ierr);
- *ynorm = snes->maxstep;
+ ierr = VecScale(Y,maxstep/(ynorm));CHKERRQ(ierr);
+ ynorm = maxstep;
}
- ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr);
+ ierr = VecMaxPointwiseDivide(Y,X,&rellength);CHKERRQ(ierr);
+#if defined(PYLITH_CUSTOM_LINESEARCH)
// Place reasonable absolute limit on minimum lambda
- minlambda = std::max(snes->steptol/rellength, 1.0/PYLITH_MAXSCALAR);
+ minlambda = std::max(steptol/rellength, 1.0/PYLITH_MAXSCALAR);
+#else // ORIGINAL
+ minlambda = steptol/rellength;
+#endif
- ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr);
+ ierr = MatMult(jac,Y,W);CHKERRQ(ierr);
#if defined(PETSC_USE_COMPLEX)
- ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr);
+ ierr = VecDot(F,W,&cinitslope);CHKERRQ(ierr);
initslope = PetscRealPart(cinitslope);
#else
- ierr = VecDot(f,w,&initslope);CHKERRQ(ierr);
+ ierr = VecDot(F,W,&initslope);CHKERRQ(ierr);
#endif
if (initslope > 0.0) initslope = -initslope;
if (initslope == 0.0) initslope = -1.0;
- ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr);
+ ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
+ if (linesearch->ops->viproject) {
+ ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
+ }
if (snes->nfuncs >= snes->max_funcs) {
ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr);
- *flag = PETSC_FALSE;
snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
- goto theend1;
- }
- ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
- if (snes->domainerror) {
- ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+ ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
- ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr);
- if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
- ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n",(double)fnorm,(double)*gnorm);CHKERRQ(ierr);
- if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + snes->ls_alpha*initslope) { /* Sufficient reduction */
- if (snes->ls_monitor) {
- ierr = PetscViewerASCIIAddTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
- ierr = PetscViewerASCIIPrintf(snes->ls_monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n",(double)fnorm,(double)*gnorm);CHKERRQ(ierr);
- ierr = PetscViewerASCIISubtractTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
- }
- goto theend1;
- }
-
- /* Fit points with quadratic */
- lambda = 1.0;
- lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope);
- lambdaprev = lambda;
- gnormprev = *gnorm;
- if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
- if (lambdatemp <= .1*lambda) lambda = .1*lambda;
- else lambda = lambdatemp;
-
- ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr);
- if (snes->nfuncs >= snes->max_funcs) {
- ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr);
- *flag = PETSC_FALSE;
- snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
- goto theend1;
- }
- ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
- if (snes->domainerror) {
- ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+ ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
+ ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
+ if (domainerror) {
+ ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
- ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr);
- if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
- if (snes->ls_monitor) {
- ierr = PetscViewerASCIIAddTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
- ierr = PetscViewerASCIIPrintf(snes->ls_monitor," Line search: gnorm after quadratic fit %14.12e\n",(double)*gnorm);CHKERRQ(ierr);
- ierr = PetscViewerASCIISubtractTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
+ if (linesearch->ops->vinorm) {
+ gnorm = fnorm;
+ ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
+ } else {
+ ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
}
- if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*snes->ls_alpha*initslope) { /* sufficient reduction */
- if (snes->ls_monitor) {
- ierr = PetscViewerASCIIAddTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
- ierr = PetscViewerASCIIPrintf(snes->ls_monitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr);
- ierr = PetscViewerASCIISubtractTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
- }
- goto theend1;
- }
- /* Fit points with cubic */
- count = 1;
- while (PETSC_TRUE) {
- if (lambda <= minlambda) {
- if (snes->ls_monitor) {
- ierr = PetscViewerASCIIAddTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
- ierr = PetscViewerASCIIPrintf(snes->ls_monitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr);
- ierr = PetscViewerASCIIPrintf(snes->ls_monitor," Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",(double)fnorm,(double)*gnorm,(double)*ynorm,(double)minlambda,(double)lambda,(double)initslope);CHKERRQ(ierr);
- ierr = PetscViewerASCIISubtractTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
- }
-#if 0 // DEBUGGING
- assert(lsctx);
- Formulation* formulation = (Formulation*) lsctx;
- assert(formulation);
- formulation->printState(&w, &g, &x, &y);
-#endif
-#if 0 // Original PETSc code
- *flag = PETSC_FALSE; // DIVERGED_LINE_SEARCH
-#else
- std::cerr << "WARNING: Line search diverged ... continuing nonlinear iterations anyway in hopes that solution will converge anyway."
- << std::endl;
-#endif
- break;
+ if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
+ ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr);
+ if (.5*gnorm*gnorm <= .5*fnorm*fnorm + lambda*alpha*initslope) { /* Sufficient reduction */
+ if (monitor) {
+ ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ ierr = PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr);
+ ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
}
- t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope;
- t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope;
- a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
- b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
- d = b*b - 3*a*initslope;
- if (d < 0.0) d = 0.0;
- if (a == 0.0) {
- lambdatemp = -initslope/(2.0*b);
- } else {
- lambdatemp = (-b + PetscSqrtScalar(d))/(3.0*a);
- } // if/else
+ } else {
+ /* Fit points with quadratic */
+ lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*lambda*initslope);
lambdaprev = lambda;
- gnormprev = *gnorm;
- if (lambdatemp > .5*lambda)
- lambdatemp = .5*lambda;
- if (lambdatemp <= .1*lambda)
- lambda = .1*lambda;
- else
- lambda = lambdatemp;
+ gnormprev = gnorm;
+ if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
+ if (lambdatemp <= .1*lambda) lambda = .1*lambda;
+ else lambda = lambdatemp;
- ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr);
+ ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
+ if (linesearch->ops->viproject) {
+ ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
+ }
if (snes->nfuncs >= snes->max_funcs) {
- ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr);
- ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr);
- *flag = PETSC_FALSE;
+ ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr);
snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
- break;
+ ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
}
- ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
- if (snes->domainerror) {
- ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+ ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
+ ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
+ if (domainerror) {
PetscFunctionReturn(0);
}
- ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr);
- if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
- if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*snes->ls_alpha*initslope) { /* is reduction enough? */
- if (snes->ls_monitor) {
- ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)*gnorm,(double)lambda);CHKERRQ(ierr);
+ if (linesearch->ops->vinorm) {
+ gnorm = fnorm;
+ ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
+ } else {
+ ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
+ }
+ if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
+ if (monitor) {
+ ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ ierr = PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",gnorm);CHKERRQ(ierr);
+ ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ }
+ if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
+ if (monitor) {
+ ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr);
+ ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
}
- break;
} else {
- if (snes->ls_monitor) {
- ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)*gnorm,(double)lambda);CHKERRQ(ierr);
+ /* Fit points with cubic */
+ for (count = 0; count < max_its; count++) {
+ if (lambda <= minlambda) {
+ if (monitor) {
+ ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ ierr = PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr);
+ ierr = PetscViewerASCIIPrintf(monitor,
+ " Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
+ fnorm, gnorm, ynorm, minlambda, lambda, initslope);CHKERRQ(ierr);
+ ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ }
+#if defined(PYLITH_CUSTOM_LINESEARCH)
+#if 0 // DEBUGGING
+ assert(lsctx);
+ Formulation* formulation = (Formulation*) lsctx;
+ assert(formulation);
+ formulation->printState(&w, &g, &x, &y);
+ std::cerr << "WARNING: Line search diverged ... continuing nonlinear iterations anyway in hopes that solution will converge anyway."
+ << std::endl;
+#endif
+ break;
+#else // ORIGINAL
+ ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+#endif
+ }
+ if (bt->order == SNES_LINESEARCH_BT_CUBIC) {
+ t1 = .5*(gnorm*gnorm - fnorm*fnorm) - lambda*initslope;
+ t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope;
+ a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
+ b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
+ d = b*b - 3*a*initslope;
+ if (d < 0.0) d = 0.0;
+ if (a == 0.0) {
+ lambdatemp = -initslope/(2.0*b);
+ } else {
+ lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
+ }
+ } else if (bt->order == SNES_LINESEARCH_BT_QUADRATIC) {
+ lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*initslope);
+ }
+ lambdaprev = lambda;
+ gnormprev = gnorm;
+ if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
+ if (lambdatemp <= .1*lambda) lambda = .1*lambda;
+ else lambda = lambdatemp;
+ ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
+ if (snes->nfuncs >= snes->max_funcs) {
+ ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr);
+ ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
+ fnorm,gnorm,ynorm,lambda,initslope);CHKERRQ(ierr);
+ ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
+ snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
+ PetscFunctionReturn(0);
+ }
+ ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
+ ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
+ if (domainerror) {
+ ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+ }
+ if (linesearch->ops->vinorm) {
+ gnorm = fnorm;
+ ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
+ } else {
+ ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
+ }
+ if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
+ if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */
+ if (monitor) {
+ ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ if (bt->order == SNES_LINESEARCH_BT_CUBIC) {
+ ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr);
+ } else {
+ ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr);
+ }
+ ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ }
+ break;
+ } else {
+ if (monitor) {
+ ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ if (bt->order == SNES_LINESEARCH_BT_CUBIC) {
+ ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr);
+ } else {
+ ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr);
+ }
+ ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
+ }
+ }
}
}
- count++;
}
- theend1:
- /* Optional user-defined check for line search step validity */
- if (snes->ops->postcheckstep && *flag) {
- ierr = (*snes->ops->postcheckstep)(snes,x,y,w,snes->postcheck,&changed_y,&changed_w);CHKERRQ(ierr);
- if (changed_y) {
- ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr);
+
+ /* postcheck */
+ ierr = SNESLineSearchPostCheck(linesearch, &changed_y, &changed_w);CHKERRQ(ierr);
+ if (changed_y) {
+ ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
+ if (linesearch->ops->viproject) {
+ ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
}
- if (changed_y || changed_w) { /* recompute the function if the step has changed */
- ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
- if (snes->domainerror) {
- ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
- PetscFunctionReturn(0);
- }
- ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr);
- if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
- ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr);
- ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr);
- ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr);
+ }
+#if defined(PYLITH_CUSTOM_LINESEARCH)
+ if (true) {
+#else // ORIGINAL
+ if (changed_y || changed_w) { /* recompute the function if the step has changed */
+#endif
+ ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
+ ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
+ if (domainerror) {
+ ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
}
- }
- ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+ if (linesearch->ops->vinorm) {
+ gnorm = fnorm;
+ ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
+ } else {
+ ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
+ }
+ ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr);
+ if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
- // ======================================================================
- ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
- if (snes->domainerror) {
- ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
- PetscFunctionReturn(0);
}
- ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr);
- if (PetscIsInfOrNanReal(*gnorm))
- SETERRQ(PETSC_COMM_SELF,
- PETSC_ERR_FP, "User provided compute function generated a "
- "Not-a-Number");
- ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr);
- ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr);
- ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr);
- // ======================================================================
+ /* copy the solution over */
+ ierr = VecCopy(W, X);CHKERRQ(ierr);
+ ierr = VecCopy(G, F);CHKERRQ(ierr);
+ ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr);
+ ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr);
+ ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr);
PetscFunctionReturn(0);
} // lineSearch
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.hh 2012-03-19 18:14:13 UTC (rev 19812)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.hh 2012-03-19 18:15:38 UTC (rev 19813)
@@ -112,34 +112,13 @@
/** Generic C interface for customized PETSc line search.
*
- * @param snes PETSc SNES solver.
+ * @param linesearch PETSc line search.
* @param lsctx Optional context for line search (not used here).
- * @param x Current iterate.
- * @param f Residual evaluated at x.
- * @param y Search direction.
- * @param w Work vector
- * @param f 2-norm of f.
- * @param xnorm Norm of x if known, otherwise 0.
- * @param g Residual evaluated at new iterate y.
- * @param w New iterate.
- * @param gnorm 2-norm of g.
- * @param ynorm 2-norm of search length.
- * @param PETSC_TRUE if line search succeeds; PETSC_FALSE on failure.
* @returns PETSc error code.
*/
static
- PetscErrorCode lineSearch(PetscSNES snes,
- void *lsctx,
- PetscVec x,
- PetscVec f,
- PetscVec y,
- PetscReal fnorm,
- PetscReal xnorm,
- PetscVec g,
- PetscVec w,
- PetscReal *ynorm,
- PetscReal *gnorm,
- PetscBool *flag);
+ PetscErrorCode lineSearch(PetscSNESLineSearch linesearch,
+ void *lsctx);
/** Generic C interface for customized PETSc initial guess.
*
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/utils/petscfwd.h
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/utils/petscfwd.h 2012-03-19 18:14:13 UTC (rev 19812)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/utils/petscfwd.h 2012-03-19 18:15:38 UTC (rev 19813)
@@ -43,6 +43,9 @@
/// forward declaration for PETSc SNES
typedef struct _p_SNES* PetscSNES;
+/// forward declatation for PETSc line search
+typedef struct _p_LineSearch* PetscSNESLineSearch;
+
/// forward declaration for PETSc PC
typedef struct _p_PC* PetscPC;
More information about the CIG-COMMITS
mailing list