C++ SolverState::finalState方法代码示例

int Epetra_TSFOperator::ApplyInverse(const Epetra_MultiVector& in, Epetra_MultiVector& out) const
{
  
  TEST_FOR_EXCEPTION(solver_.ptr().get()==0, std::runtime_error,
		     "no solver provided for Epetra_TSFOperator::ApplyInverse");
  TEST_FOR_EXCEPTION(!isNativeEpetra_ && !isCompoundEpetra_, std::runtime_error,
		     "Epetra_TSFOperator::ApplyInverse expects either "
		     "a native epetra operator or a compound operator with "
		     "Epetra domain and range spaces");
  const Epetra_Vector* cevIn = dynamic_cast<const Epetra_Vector*>(&in);
  Epetra_Vector* evIn = const_cast<Epetra_Vector*>(cevIn);
  Epetra_Vector* evOut = dynamic_cast<Epetra_Vector*>(&out);

  TEST_FOR_EXCEPTION(evIn==0, std::runtime_error, "Epetra_TSFOperator::Apply "
		     "cannot deal with multivectors");
  TEST_FOR_EXCEPTION(evOut==0, std::runtime_error, "Epetra_TSFOperator::Apply "
		     "cannot deal with multivectors");

  const EpetraVectorSpace* ed 
    = dynamic_cast<const EpetraVectorSpace*>(A_.range().ptr().get());
  const EpetraVectorSpace* er 
    = dynamic_cast<const EpetraVectorSpace*>(A_.domain().ptr().get());

  RCP<Thyra::VectorBase<double> > vpIn 
    = rcp(new EpetraVector(rcp(ed, false), 
			   rcp(evIn, false)));
  RCP<Thyra::VectorBase<double> > vpOut 
    = rcp(new EpetraVector(rcp(er, false), 
			   rcp(evOut, false)));
  Vector<double> vIn = vpIn;
  Vector<double> vOut = vpOut;
  
  SolverState<double> state = solver_.solve(A_, vIn, vOut);

  if (state.finalState() == SolveCrashed) return -1;
  else if (state.finalState() == SolveFailedToConverge) return -2;
  else out = EpetraVector::getConcrete(vOut);

  return 0;
}

void LinearPDEConstrainedObj::solveState(const Vector<double>& x) const
{
  Tabs tab(0);
  PLAYA_MSG2(verb(), tab << "solving state"); 
  PLAYA_MSG3(verb(), tab << "|x|=" << x.norm2()); 
  PLAYA_MSG5(verb(), tab << "x=" << endl << tab << x.norm2());
  setDiscreteFunctionVector(designVarVal(), x);

  /* solve the state equations in order */
  for (int i=0; i<stateProbs_.size(); i++)
  {
    SolverState<double> status 
      = stateProbs_[i].solve(solvers_[i], stateVarVals(i));
    TEUCHOS_TEST_FOR_EXCEPTION(status.finalState() != SolveConverged,
      std::runtime_error,
      "state equation could not be solved: status="
      << status.stateDescription());
  }

  PLAYA_MSG2(verb(), tab << "done state solve"); 
  /* do postprocessing */
  statePostprocCallback();
}

  template <class Scalar> inline
  SolverState<Scalar> BlockTriangularSolver<Scalar>
  ::solve(const LinearOperator<Scalar>& op,
          const Vector<Scalar>& rhs,
          Vector<Scalar>& soln) const
  {
    int nRows = op.numBlockRows();
    int nCols = op.numBlockCols();

    soln = op.domain().createMember();
    //    bool converged = false;

    TEST_FOR_EXCEPTION(nRows != rhs.space().numBlocks(), std::runtime_error,
                       "number of rows in operator " << op
                       << " not equal to number of blocks on RHS "
                       << rhs);

    TEST_FOR_EXCEPTION(nRows != nCols, std::runtime_error,
                       "nonsquare block structure in block triangular "
                       "solver: nRows=" << nRows << " nCols=" << nCols);

    bool isUpper = false;
    bool isLower = false;

    for (int r=0; r<nRows; r++)
      {
        for (int c=0; c<nCols; c++)
          {
            if (op.getBlock(r,c).ptr().get() == 0 ||
                dynamic_cast<const SimpleZeroOp<Scalar>* >(op.getBlock(r,c).ptr().get()))
              {
                TEST_FOR_EXCEPTION(r==c, std::runtime_error,
                                   "zero diagonal block (" << r << ", " << c 
                                   << " detected in block "
                                   "triangular solver. Operator is " << op);
                continue;
              }
            else
              {
                if (r < c) isUpper = true;
                if (c < r) isLower = true;
              }
          }
      }

    TEST_FOR_EXCEPTION(isUpper && isLower, std::runtime_error, 
                       "block triangular solver detected non-triangular operator "
                       << op);

    bool oneSolverFitsAll = false;
    if ((int) solvers_.size() == 1 && nRows != 1) 
      {
        oneSolverFitsAll = true;
      }

    for (int i=0; i<nRows; i++)
      {
        int r = i;
        if (isUpper) r = nRows - 1 - i;
        Vector<Scalar> rhs_r = rhs.getBlock(r);
        for (int j=0; j<i; j++)
          {
            int c = j;
            if (isUpper) c = nCols - 1 - j;
            if (op.getBlock(r,c).ptr().get() != 0)
              {
                rhs_r = rhs_r - op.getBlock(r,c) * soln.getBlock(c);
              }
          }

        SolverState<Scalar> state;
        Vector<Scalar> soln_r;
        if (oneSolverFitsAll)
          {
            state = solvers_[0].solve(op.getBlock(r,r), rhs_r, soln_r);
          }
        else
          {
            state = solvers_[r].solve(op.getBlock(r,r), rhs_r, soln_r);
          }
        if (nRows > 1) soln.setBlock(r, soln_r);
        else soln = soln_r;
        if (state.finalState() != SolveConverged)
          {
            return state;
          }
      }

    return SolverState<Scalar>(SolveConverged, "block solves converged",
                               0, ScalarTraits<Scalar>::zero());
  }

int main(int argc, char *argv[]) 
{
  try
  {
    GlobalMPISession session(&argc, &argv);


    MPIComm::world().synchronize();

    /* create the nonlinear operator */
    VectorType<double> type = new EpetraVectorType();
    int nProc = MPIComm::world().getNProc();
    int nLocalRows = 128/nProc;
    PoissonBoltzmannOp* prob = new PoissonBoltzmannOp(nLocalRows, type);
    NonlinearOperator<double> F = prob;

    /* create the nox solver */
    ParameterXMLFileReader reader("nox.xml");
    ParameterList noxParams = reader.getParameters();

    Out::root() << "solver params = " << noxParams << std::endl;

    NOXSolver solver(noxParams);

    Vector<double> soln;
    SolverState<double> stat = solver.solve(F, soln);
    TEUCHOS_TEST_FOR_EXCEPTION(stat.finalState() != SolveConverged,
      runtime_error, "solve failed");

    Out::root() << "numerical solution = " << std::endl;
    Out::os() << soln << std::endl;

    Vector<double> exact = prob->exactSoln();

    Out::root() << "exact solution = " << std::endl;
    Out::os() << exact << std::endl;

//bvbw reddish port hack
    double temp_val = nLocalRows*nProc;
    double err = (exact-soln).norm2()/sqrt(temp_val);
    Out::root() << "error norm = " << err << std::endl;
      

    double tol = 1.0e-6;
    if (err > tol)
    {
      Out::root() << "NOX Poisson-Boltzmann test FAILED" << std::endl;
      return 1;
    }
    else
    {
      Out::root() << "NOX Poisson-Boltzmann test PASSED" << std::endl;
      return 0;
    }
  }
  catch(std::exception& e)
  {
    Out::root() << "Caught exception: " << e.what() << std::endl;
    return -1;
  }
}

int main(int argc, char** argv)
{
  
  try
		{
      Sundance::init(&argc, &argv);
            
      /* We will do our linear algebra using Epetra */
      VectorType<double> vecType = new EpetraVectorType();

      /* Create a mesh. It will be of type BasisSimplicialMesh, and will
       * be built using a PartitionedLineMesher. */
      MeshType meshType = new BasicSimplicialMeshType();
      int nx = 32;

      MeshSource mesher = new PartitionedRectangleMesher(
        0.0, 1.0, nx,
        0.0, 1.0, nx,
        meshType);
      Mesh mesh = mesher.getMesh();

      /* Create a cell filter that will identify the maximal cells
       * in the interior of the domain */
      CellFilter interior = new MaximalCellFilter();

      /* Make cell filters for the east and west boundaries */
      CellFilter edges = new DimensionalCellFilter(1);
      CellFilter west = edges.coordSubset(0, 0.0);
      CellFilter east = edges.coordSubset(0, 1.0);

      /* Create unknown and test functions */
      Expr u = new UnknownFunction(new Lagrange(1), "u");
      Expr v = new TestFunction(new Lagrange(1), "v");

      /* Create differential operator and coordinate function */
      Expr x = new CoordExpr(0);
      Expr grad = gradient(1);

      /* We need a quadrature rule for doing the integrations */
      QuadratureFamily quad = new GaussianQuadrature(4);

      /* Define the parameter */
      Expr xi = new Sundance::Parameter(0.0);

      /* Construct a forcing term to provide an exact solution for
       * validation purposes. This will involve the parameter. */
      Expr uEx = x*(1.0-x)*(1.0+xi*exp(x));
      Expr f = -(-20 - exp(x)*xi*(1 + 32*x + 10*x*x + 
          exp(x)*(-1 + 2*x*(2 + x))*xi))/10.0;

      /* Define the weak form, using the parameter expression. This weak form
       * can be used for all parameter values. */
      Expr eqn = Integral(interior, 
        (1.0 + 0.1*xi*exp(x))*(grad*v)*(grad*u) - f*v, quad);

      /* Define the Dirichlet BC */
      Expr h = new CellDiameterExpr();
      Expr bc = EssentialBC(east + west, v*u/h, quad);

      /* We can now set up the linear problem. This can be reused
       * for different parameter values. */
      LinearProblem prob(mesh, eqn, bc, v, u, vecType);

      /* make a projector for the exact solution. Just like the
       * problem, this can be reused for different parameter values. */
      DiscreteSpace ds(mesh, new Lagrange(1), vecType);
      L2Projector proj(ds, uEx);

      /* Get the solver and declare variables for the results */
      LinearSolver<double> solver = LinearSolverBuilder::createSolver("aztec-ml.xml");
      Expr soln;
      SolverState<double> state;

      /* Set up the sweep from xi=0 to xi=xiMax in nSteps steps. */
      int nSteps = 10;
      double xiMax = 2.0;
      
      /* Make an array in which to keep the observed errors */
      Array<double> err(nSteps);

      /* Do the sweep */
      for (int n=0; n<nSteps; n++)
      {
        /* Update the parameter value */
        double xiVal = xiMax*n/(nSteps - 1.0);
        xi.setParameterValue(xiVal);
        Out::root() << "step n=" << n << " of " << nSteps << " xi=" << xiVal;

        /* Solve the problem. The updated parameter value is automatically used. */
        state = prob.solve(solver, soln);

        TEUCHOS_TEST_FOR_EXCEPTION(state.finalState() != SolveConverged,
          std::runtime_error,
          "solve failed!");

        /* Project the exact solution onto a discrrete space for viz. The updated
         * parameter value is automatically used. */
        Expr uEx0 = proj.project();

        /* Write the approximate and exact solutions for viz */
//.........这里部分代码省略.........

void LinearPDEConstrainedObj
::solveStateAndAdjoint(const Vector<double>& x) const
{
  Tabs tab(0);
  PLAYA_MSG2(verb(), tab << "solving state and adjoint"); 
  PLAYA_MSG3(verb(), tab << "|x|=" << x.norm2()); 
  PLAYA_MSG5(verb(), tab << "x=" << endl << tab << x.norm2()); 

  Tabs tab1;
  setDiscreteFunctionVector(designVarVal(), x);

  PLAYA_MSG3(verb(), tab1 << "solving state eqns");
  /* solve the state equations in order */
  for (int i=0; i<stateProbs_.size(); i++)
  {
    SolverState<double> status 
      = stateProbs_[i].solve(solvers_[i], stateVarVals(i));

    /* if the solve failed, write out the design var and known state
     * variables */
    if (status.finalState() != SolveConverged)
    {
      FieldWriter w = new VTKWriter("badSolve");
      w.addMesh(Lagrangian().mesh());
      w.addField("designVar", new ExprFieldWrapper(designVarVal()));
      for (int j=0; j<i; j++)
      {
        Expr tmp = stateVarVals(j).flatten();
        for (int k=0; k<tmp.size(); k++)
        {
          w.addField("stateVar-"+Teuchos::toString(j)+"-"+Teuchos::toString(k),
            new ExprFieldWrapper(tmp[k]));
        }
      }
      w.write();
    }
    TEUCHOS_TEST_FOR_EXCEPTION(status.finalState() != SolveConverged,
      std::runtime_error,
      "state equation " << i 
      << " could not be solved: status="
      << status.stateDescription());
  }

  PLAYA_MSG3(verb(), tab1 << "done solving state eqns");

  /* do postprocessing */
  statePostprocCallback();

  PLAYA_MSG3(verb(), tab1 << "solving adjoint eqns");

  /* solve the adjoint equations in reverse order */
  for (int i=adjointProbs_.size()-1; i>=0; i--)
  {
    SolverState<double> status 
      = adjointProbs_[i].solve(solvers_[i], adjointVarVals(i));

    /* if the solve failed, write out the design var and known state
     * and adjoint variables */
    if (status.finalState() != SolveConverged)
    {
      FieldWriter w = new VTKWriter("badSolve");
      w.addMesh(Lagrangian().mesh());
      w.addField("designVar", new ExprFieldWrapper(designVarVal()));
      for (int j=0; j<stateProbs_.size(); j++)
      {
        Expr tmp = stateVarVals(j).flatten();
        for (int k=0; k<tmp.size(); k++)
        {
          w.addField("stateVar-"+Teuchos::toString(j)+"-"+Teuchos::toString(k),
            new ExprFieldWrapper(tmp[k]));
        }
      }
      for (int j=adjointProbs_.size()-1; j>i; j--)
      {
        Expr tmp = adjointVarVals(j).flatten();
        for (int k=0; k<tmp.size(); k++)
        {
          w.addField("adjointVar-"+Teuchos::toString(j)+"-"+Teuchos::toString(k),
            new ExprFieldWrapper(tmp[k]));
        }

      }
      w.write();

    }
    TEUCHOS_TEST_FOR_EXCEPTION(status.finalState() != SolveConverged,
      std::runtime_error,
      "adjoint equation " << i 
      << " could not be solved: status="
      << status.stateDescription());
  }
  PLAYA_MSG3(verb(), tab1 << "done solving adjoint eqns");
  PLAYA_MSG2(verb(), tab1 << "done solving state and adjoint eqns");
}

int main(int argc, char** argv)
{
  try
  {
    int nx = 32;
    double convTol = 1.0e-8;
    double lambda = 0.5;
    Sundance::setOption("nx", nx, "Number of elements");
    Sundance::setOption("tol", convTol, "Convergence tolerance");
    Sundance::setOption("lambda", lambda, "Lambda (parameter in Bratu's equation)");

    Sundance::init(&argc, &argv);

    Out::root() << "Bratu problem (lambda=" << lambda << ")" << endl;
    Out::root() << "Newton's method with automated linearization" 
                << endl << endl;

    VectorType<double> vecType = new EpetraVectorType();

    MeshType meshType = new BasicSimplicialMeshType();
    MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, nx, meshType);
    Mesh mesh = mesher.getMesh();

    CellFilter interior = new MaximalCellFilter();
    CellFilter sides = new DimensionalCellFilter(mesh.spatialDim()-1);
    CellFilter left = sides.subset(new CoordinateValueCellPredicate(0, 0.0));
    CellFilter right = sides.subset(new CoordinateValueCellPredicate(0, 1.0));
    
    BasisFamily basis = new Lagrange(1);
    Expr u = new UnknownFunction(basis, "w");
    Expr v = new TestFunction(basis, "v");

    Expr grad = gradient(1);

    Expr x = new CoordExpr(0);

    const double pi = 4.0*atan(1.0);
    Expr uExact = sin(pi*x);
    Expr R = pi*pi*uExact - lambda*exp(uExact);

    QuadratureFamily quad4 = new GaussianQuadrature(4);
    QuadratureFamily quad2 = new GaussianQuadrature(2);

    DiscreteSpace discSpace(mesh, basis, vecType);
    Expr uPrev = new DiscreteFunction(discSpace, 0.5);

    Expr eqn 
      = Integral(interior, (grad*v)*(grad*u) - v*lambda*exp(u) - v*R, quad4);

    Expr h = new CellDiameterExpr();
    Expr bc = EssentialBC(left+right, v*u/h, quad2); 

    NonlinearProblem prob(mesh, eqn, bc, v, u, uPrev, vecType);

    NonlinearSolver<double> solver 
      = NonlinearSolverBuilder::createSolver("playa-newton-amesos.xml");

    Out::root() << "Newton solve" << endl;

    SolverState<double> state = prob.solve(solver);
    
    TEUCHOS_TEST_FOR_EXCEPTION(state.finalState() != SolveConverged,
      std::runtime_error,
      "Nonlinear solve failed to converge: message=" << state.finalMsg());
    
    Expr soln = uPrev;
    FieldWriter writer = new DSVWriter("AutoLinearizedBratu.dat");
    writer.addMesh(mesh);
    writer.addField("soln", new ExprFieldWrapper(soln[0]));
    writer.write();

    Out::root() << "Converged!" << endl << endl;

    double L2Err = L2Norm(mesh, interior, soln-uExact, quad4);
    Out::root() << "L2 Norm of error: " << L2Err << endl;
    
    Sundance::passFailTest(L2Err, 1.5/((double) nx*nx));
  }
	catch(std::exception& e) 
  {
    Sundance::handleException(e);
  }
  Sundance::finalize(); 
  return Sundance::testStatus();
}