<p><b>cnewman@lanl.gov</b> 2010-08-06 15:51:07 -0600 (Fri, 06 Aug 2010)</p><p>More additions to implicit branch<br>
</p><hr noshade><pre><font color="gray">Modified: branches/implicit/Makefile
===================================================================
--- branches/implicit/Makefile 2010-08-06 19:14:48 UTC (rev 469)
+++ branches/implicit/Makefile 2010-08-06 21:51:07 UTC (rev 470)
@@ -159,6 +159,7 @@
#########################
ifdef TRILINOS_HOME
include $(TRILINOS_HOME)/include/Makefile.export.Epetra
+ include $(TRILINOS_HOME)/include/Makefile.export.NOX
MPI_HOME=/opt/OpenMPI/openmpi-1.2.8-gcc/ib
MPIINC = -I$(MPI_HOME)/include
MPILIB = -L$(MPI_HOME)/lib64 -lmpi_cxx
@@ -173,7 +174,7 @@
ifdef TRILINOS_LIB_PATH
LIBS += -L$(TRILINOS_LIB_PATH) -lzoltan
else
- LIBS += -L$(TRILINOS_HOME)/lib $(EPETRA_LIBRARIES) $(MPILIB) $(EPETRA_TPL_LIBRARIES) -lstdc++
+ LIBS += -L$(TRILINOS_HOME)/lib $(EPETRA_LIBRARIES) $(MPILIB) $(EPETRA_TPL_LIBRARIES) $(NOX_LIBRARIES) -lstdc++
endif
endif
#########################
Modified: branches/implicit/src/core_sw/Makefile
===================================================================
--- branches/implicit/src/core_sw/Makefile 2010-08-06 19:14:48 UTC (rev 469)
+++ branches/implicit/src/core_sw/Makefile 2010-08-06 21:51:07 UTC (rev 470)
@@ -6,7 +6,7 @@
module_time_integration_support.o \
module_implicit_integration.o \
module_time_integration.o \
- mpas_interface.o nox_driver.o
+ mpas_interface.o nox_driver.o nox_interface.o
all: core_sw
@@ -19,7 +19,7 @@
module_implicit_integration.o:
-mpas_interface.o: module_time_integration_support.o module_implicit_integration.o module_test_cases.o module_time_integration.o nox_driver.o
+mpas_interface.o: module_time_integration_support.o module_implicit_integration.o module_test_cases.o module_time_integration.o nox_driver.o nox_interface.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/implicit/src/core_sw/nox_driver.C
===================================================================
--- branches/implicit/src/core_sw/nox_driver.C 2010-08-06 19:14:48 UTC (rev 469)
+++ branches/implicit/src/core_sw/nox_driver.C 2010-08-06 21:51:07 UTC (rev 470)
@@ -1,6 +1,6 @@
-
#include "NOX.H"
#include "NOX_Epetra.H"
+#include "Epetra_Vector.h"
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
@@ -31,6 +31,9 @@
}
static Teuchos::RCP<Problem_Interface> interface;
+static Teuchos::RCP<NOX::Epetra::LinearSystem> linsys;
+static Teuchos::RCP<NOX::Solver::Generic> solver;
+static Teuchos::RCP<NOX::Epetra::Vector> initialGuess;
extern "C" {
@@ -120,8 +123,77 @@
// Create the interface between the test problem and the nonlinear solver
// This is created by the user using inheritance of the abstract base class:
// NOX_Epetra_Interface
- //interface = Teuchos::rcp(new Problem_Interface(*nelems, statevector, Comm,
- // blackbox_res, residualFunction));
+ interface = Teuchos::rcp(new Problem_Interface(*nelems, statevector, Comm,
+ blackbox_res, residualFunction));
+
+ Teuchos::RCP<NOX::Epetra::Interface::Required> iReq = interface;
+
+ // Create the Epetra_Vector for the state vector
+ //Teuchos::RCP<Epetra_Vector> soln=Teuchos::rcp(new Epetra_Vector(interface->getVector()));
+ Teuchos::RCP<Epetra_Vector> soln = interface->getVector();
+
+ // Create the nox vector
+ initialGuess = Teuchos::RCP<NOX::Epetra::Vector>(new NOX::Epetra::Vector(soln));
+
+ Teuchos::RCP<NOX::Epetra::MatrixFree> FD =
+ Teuchos::rcp(new NOX::Epetra::MatrixFree(nlPrintParams, interface, soln));
+
+ Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac = FD;
+
+// note: it doesn't matter if this is null
+ Teuchos::RCP<NOX::Epetra::Scaling> scaling = interface->getScaling();
+// Create the Group
+ Teuchos::RCP<NOX::Epetra::Group> grp =
+ Teuchos::rcp(new NOX::Epetra::Group(nlPrintParams, iReq, *initialGuess,
+ linsys));
+
+// Create the convergence tests
+ Teuchos::RCP<NOX::StatusTest::NormF> absresid =
+ Teuchos::rcp(new NOX::StatusTest::NormF(1.0e-8));
+ Teuchos::RCP<NOX::StatusTest::NormF> relresid =
+ Teuchos::rcp(new NOX::StatusTest::NormF(*grp.get(), 1.0e-4));
+ Teuchos::RCP<NOX::StatusTest::NormUpdate> update =
+ Teuchos::rcp(new NOX::StatusTest::NormUpdate(1.0e-5));
+ Teuchos::RCP<NOX::StatusTest::NormWRMS> wrms =
+ Teuchos::rcp(new NOX::StatusTest::NormWRMS(1.0e-2, 1.0e-8));
+ Teuchos::RCP<NOX::StatusTest::Combo> converged =
+ Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::AND));
+ converged->addStatusTest(absresid);
+ converged->addStatusTest(relresid);
+ converged->addStatusTest(wrms);
+ converged->addStatusTest(update);
+ Teuchos::RCP<NOX::StatusTest::MaxIters> maxiters =
+ Teuchos::rcp(new NOX::StatusTest::MaxIters(200)); // cut this back for now
+// Teuchos::rcp(new NOX::StatusTest::MaxIters(20));
+ Teuchos::RCP<NOX::StatusTest::FiniteValue> fv =
+ Teuchos::rcp(new NOX::StatusTest::FiniteValue);
+ Teuchos::RCP<NOX::StatusTest::Combo> combo =
+ Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
+ combo->addStatusTest(fv);
+ combo->addStatusTest(converged);
+ combo->addStatusTest(maxiters);
+ Teuchos::RCP<Teuchos::ParameterList> nlParamsRCP = Teuchos::rcp(&nlParams, false);
+ solver = NOX::Solver::buildSolver(grp, combo, nlParamsRCP);
+
}
+void noxsolve(int* nelems, double* statevector,
+ void* blackbox_res, void* blackbox_prec)
+{
+ // reset solver with new initial Guess
+ Epetra_Vector& initialGuessEV = initialGuess->getEpetraVector();
+ for (int i=0; i<*nelems; i++) initialGuessEV[i] = statevector[i];
+ solver->reset(*initialGuess);
+ // reset interface with new blackbox
+ interface->resetBlackbox(blackbox_res, blackbox_prec);
+ NOX::StatusTest::StatusType status = solver->solve();
+
+ // copy out final solution
+ const Epetra_Vector& finalSolution = (dynamic_cast<const NOX::Epetra::Vector&>
+ (solver->getSolutionGroup().getX())).getEpetraVector();
+ for (int i=0; i<*nelems; i++) statevector[i] = finalSolution[i];
+
}
+
+
+}
Modified: branches/implicit/src/core_sw/nox_interface.C
===================================================================
--- branches/implicit/src/core_sw/nox_interface.C 2010-08-06 19:14:48 UTC (rev 469)
+++ branches/implicit/src/core_sw/nox_interface.C 2010-08-06 21:51:07 UTC (rev 470)
@@ -1,15 +1,24 @@
-#include nox_interface.h
+#include "nox_interface.h"
+//using Teuchos::RCP;
+//using Teuchos::rcp;
+Problem_Interface::Problem_Interface(int nelems,
+ double* statevector,
+ const Epetra_Comm& comm_,
+ void* blackbox_res,
+ void (*residualFunction)(double *, double *, int, void *, double))
+ : N(nelems),
+ comm(comm_)
+{
+this->Create(statevector);
+}
-Problem_Interface::Problem_Interface(*nelems, statevector, Comm,
- blackbox_res, residualFunction)
+Problem_Interface::~Problem_Interface()
{ }
- ~Problem_Interface()
-{ }
//! Compute and return F
bool Problem_Interface::computeF(const Epetra_Vector& x, Epetra_Vector& FVec,
- FillType flag = Residual)
+ FillType flag)
{ }
//! Compute an explicit Jacobian
@@ -22,6 +31,45 @@
//! Computes a user supplied preconditioner based on input vector x. Returns true if computation was successful.
bool Problem_Interface::computePreconditioner(const Epetra_Vector& x, Epetra_Operator& M,
- Teuchos::ParameterList* precParams = 0)
+ Teuchos::ParameterList* precParams)
{ }
+Teuchos::RCP<Epetra_Vector> Problem_Interface::getVector() const
+{
+ return solution;
+}
+
+Teuchos::RCP<NOX::Epetra::Scaling> Problem_Interface::getScaling() const
+ {
+ return scaling;
+ }
+
+void Problem_Interface::Create(double* statevector)
+ {
+ globalMap = Teuchos::rcp(new Epetra_Map(-1, N, indexBase_, comm));
+ solution = Teuchos::rcp(new Epetra_Vector(Copy, *globalMap, statevector));
+
+// left_scaling = rcp(new Epetra_Vector(*solution));
+// right_scaling = rcp(new Epetra_Vector(*solution));
+
+// left_scaling->PutScalar(1.0);
+// right_scaling->PutScalar(1.0);
+
+// left_scaling->SetLabel("Left Scaling");
+// right_scaling->SetLabel("Right Scaling");
+
+// scaling = rcp(new NOX::Epetra::Scaling() );
+
+
+// // note: the scaling will be recomputed, but Trilinos is
+// // aware of that because it has a pointer:
+// scaling->addUserScaling(NOX::Epetra::Scaling::Left,left_scaling);
+// scaling->addUserScaling(NOX::Epetra::Scaling::Right,right_scaling);
+
+ scaling = Teuchos::null; // we don't use scaling at the moment. If we need it,
+ // we probably have to put it in the IMPOP_BlockPreconditioner
+ // as the outer iteration is Jacobian-free
+
+// label = "IMPOP Interface ("+prec_type+" preconditioning)";
+// DEBUG("Created Interface labelled '"<<label<<"'");
+ }
Modified: branches/implicit/src/core_sw/nox_interface.h
===================================================================
--- branches/implicit/src/core_sw/nox_interface.h 2010-08-06 19:14:48 UTC (rev 469)
+++ branches/implicit/src/core_sw/nox_interface.h 2010-08-06 21:51:07 UTC (rev 470)
@@ -8,20 +8,24 @@
// ---------- Standard Includes ----------
#include <iostream>
#include "Epetra_Vector.h"
+#include "Epetra_Map.h"
#include "Epetra_Operator.h"
#include "Epetra_RowMatrix.h"
#include "NOX_Epetra_Interface_Required.H" // base class
#include "NOX_Epetra_Interface_Jacobian.H" // base class
#include "NOX_Epetra_Interface_Preconditioner.H" // base class
+#include "NOX_Epetra_Scaling.H"
class Problem_Interface : public NOX::Epetra::Interface::Required,
public NOX::Epetra::Interface::Jacobian,
public NOX::Epetra::Interface::Preconditioner
{
public:
- Problem_Interface(int nelems, double* statevector,
+ Problem_Interface(int nelems,
+ double* statevector,
const Epetra_Comm& comm_,
- void* blackbox_res, void (*residualFunction)(double *, double *, int, void *, double));
+ void* blackbox_res,
+ void (*residualFunction)(double *, double *, int, void *, double));
~Problem_Interface();
//! Compute and return F
@@ -38,6 +42,27 @@
bool computePreconditioner(const Epetra_Vector& x, Epetra_Operator& M,
Teuchos::ParameterList* precParams = 0);
+ Teuchos::RCP<Epetra_Vector> getVector() const;
+ Teuchos::RCP<NOX::Epetra::Scaling> getScaling() const;
+ void resetBlackbox(void* blackbox_res_, void* blackbox_prec_) {
+ blackbox_res=blackbox_res_; blackbox_prec=blackbox_prec_; }
+
+private:
+ void* blackbox_res;
+ void* blackbox_prec;
+ Teuchos::RCP<Epetra_Vector> solution;
+ Teuchos::RCP<NOX::Epetra::Scaling> scaling;
+ //! number of elements on this process
+ int N;
+ // global communicator
+ const Epetra_Comm& comm;
+ //! index base of global map
+ static const int indexBase_=0;
+ Teuchos::RCP<Epetra_Map> globalMap;
+ //! function called by all constructors (creates map etc.)
+ void Create(double* statevector);
+
+
};
#endif
</font>
</pre>