C++ InArgs::get_p_sg方法代码示例

void
twoD_diffusion_ME::
evalModel(const InArgs& inArgs, const OutArgs& outArgs) const
{

  //
  // Determinisic calculation
  //

  // Solution vector
  Teuchos::RCP<const Epetra_Vector> det_x = inArgs.get_x();

  // Parameters
  Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(0);
  if (p == Teuchos::null)
    p = p_init;

  Teuchos::RCP<Epetra_Vector> f = outArgs.get_f();
  Teuchos::RCP<Epetra_Operator> W = outArgs.get_W();
  Teuchos::RCP<Epetra_Operator> WPrec = outArgs.get_WPrec();
  if (f != Teuchos::null || W != Teuchos::null || WPrec != Teuchos::null) {
    if (basis != Teuchos::null) {
      for (int i=0; i<point.size(); i++)
        point[i] = (*p)[i];
      basis->evaluateBases(point, basis_vals);
      A->PutScalar(0.0);
      for (int k=0;k<A_k.size();k++)
        EpetraExt::MatrixMatrix::Add((*A_k[k]), false, basis_vals[k], *A, 1.0);
    }
    else {
      *A = *(A_k[0]);
      for (int k=1;k<A_k.size();k++)
        EpetraExt::MatrixMatrix::Add((*A_k[k]), false, (*p)[k-1], *A, 1.0);
    }
    A->FillComplete();
    A->OptimizeStorage();
  }

  // Residual
  if (f != Teuchos::null) {
    Teuchos::RCP<Epetra_Vector> kx = Teuchos::rcp(new Epetra_Vector(*x_map));
    A->Apply(*det_x,*kx);
    f->Update(1.0,*kx,-1.0, *b, 0.0);
  }

  // Jacobian
  if (W != Teuchos::null) {
    Teuchos::RCP<Epetra_CrsMatrix> jac =
      Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W, true);
    *jac = *A;
    jac->FillComplete();
    jac->OptimizeStorage();
  }

  // Preconditioner
  if (WPrec != Teuchos::null)
    precFactory->recompute(A, WPrec);

  // Responses (mean value)
  Teuchos::RCP<Epetra_Vector> g = outArgs.get_g(0);
  if (g != Teuchos::null) {
    (det_x->MeanValue(&(*g)[0]));
    (*g)[0] *= double(det_x->GlobalLength()) / double(mesh.size());
  }

  //
  // Stochastic Galerkin calculation
  //

  // Stochastic solution vector
  InArgs::sg_const_vector_t x_sg = inArgs.get_x_sg();

  // Stochastic parameters
  InArgs::sg_const_vector_t p_sg = inArgs.get_p_sg(0);

  // Stochastic residual
  OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
  if (f_sg != Teuchos::null) {

    // Get stochastic expansion data
    Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expn =
      inArgs.get_sg_expansion();
    typedef Stokhos::Sparse3Tensor<int,double> Cijk_type;
    Teuchos::RCP<const Cijk_type> Cijk = expn->getTripleProduct();
    const Teuchos::Array<double>& norms = basis->norm_squared();

    if (sg_kx_vec_all.size() != basis->size()) {
      sg_kx_vec_all.resize(basis->size());
      for (int i=0;i<basis->size();i++) {
        sg_kx_vec_all[i] = Teuchos::rcp(new Epetra_Vector(*x_map));
      }
    }
    f_sg->init(0.0);

    Cijk_type::k_iterator k_begin = Cijk->k_begin();
    Cijk_type::k_iterator k_end = Cijk->k_end();
    for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
      int k = Stokhos::index(k_it);
      for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
           j_it != Cijk->j_end(k_it); ++j_it) {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                                     sacado_param_vec, p_vec.get(),
                                     NULL, NULL, NULL, NULL, g_out.get(), NULL,
                                     dgdp_out.get());
        g_computed = true;
      }
    }

    // Need to handle dg/dp for distributed p

    if (g_out != Teuchos::null && !g_computed)
      app->evaluateResponse(i, curr_time, x_dot.get(), x_dotdot.get(), *x, sacado_param_vec,
                            *g_out);
  }

  //
  // Stochastic Galerkin
  //
#ifdef ALBANY_SG_MP
  InArgs::sg_const_vector_t x_sg = inArgs.get_x_sg();
  if (x_sg != Teuchos::null) {
    app->init_sg(inArgs.get_sg_basis(),
                 inArgs.get_sg_quadrature(),
                 inArgs.get_sg_expansion(),
                 x_sg->productComm());
    InArgs::sg_const_vector_t x_dot_sg  = inArgs.get_x_dot_sg();
    InArgs::sg_const_vector_t x_dotdot_sg  = inArgs.get_x_dotdot_sg();
    if (x_dot_sg != Teuchos::null || x_dotdot_sg != Teuchos::null) {
      alpha = inArgs.get_alpha();
      omega = inArgs.get_omega();
      beta = inArgs.get_beta();
      curr_time  = inArgs.get_t();
    }

    InArgs::sg_const_vector_t epetra_p_sg = inArgs.get_p_sg(0);
    Teuchos::Array<int> p_sg_index;
    for (int i=0; i<num_param_vecs; i++) {
      InArgs::sg_const_vector_t p_sg = inArgs.get_p_sg(i);
      if (p_sg != Teuchos::null) {
        p_sg_index.push_back(i);
        for (int j=0; j<p_sg_vals[i].size(); j++) {
          int num_sg_blocks = p_sg->size();
          p_sg_vals[i][j].reset(app->getStochasticExpansion(), num_sg_blocks);
          p_sg_vals[i][j].copyForWrite();
          for (int l=0; l<num_sg_blocks; l++) {
            p_sg_vals[i][j].fastAccessCoeff(l) = (*p_sg)[l][j];
          }
        }
      }
    }

    OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
    OutArgs::sg_operator_t W_sg = outArgs.get_W_sg();
    bool f_sg_computed = false;

    // W_sg
    if (W_sg != Teuchos::null) {
      Stokhos::VectorOrthogPoly<Epetra_CrsMatrix> W_sg_crs(W_sg->basis(),
                                                           W_sg->map());
      for (int i=0; i<W_sg->size(); i++)
        W_sg_crs.setCoeffPtr(
          i,
          Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_sg->getCoeffPtr(i)));
      app->computeGlobalSGJacobian(alpha, beta, omega, curr_time,
                                   x_dot_sg.get(),  x_dotdot_sg.get(), *x_sg,
                                   sacado_param_vec, p_sg_index, p_sg_vals,
                                   f_sg.get(), W_sg_crs);

void 
Stokhos::SGQuadModelEvaluator::
evalModel(const InArgs& inArgs, const OutArgs& outArgs) const
{
  // Create underlying inargs
  InArgs me_inargs = me->createInArgs();
  if (me_inargs.supports(IN_ARG_x))
    me_inargs.set_x(inArgs.get_x());
  if (me_inargs.supports(IN_ARG_x_dot))
    me_inargs.set_x_dot(inArgs.get_x_dot());
  if (me_inargs.supports(IN_ARG_alpha))
    me_inargs.set_alpha(inArgs.get_alpha());
  if (me_inargs.supports(IN_ARG_beta))
    me_inargs.set_beta(inArgs.get_beta());
  if (me_inargs.supports(IN_ARG_t))
    me_inargs.set_t(inArgs.get_t());
  for (int i=0; i<num_p; i++)
    me_inargs.set_p(i, inArgs.get_p(i));

  // Create underlying outargs
  OutArgs me_outargs = me->createOutArgs();
  if (me_outargs.supports(OUT_ARG_f))
    me_outargs.set_f(outArgs.get_f());
  if (me_outargs.supports(OUT_ARG_W))
    me_outargs.set_W(outArgs.get_W());
  for (int j=0; j<num_p; j++)
    if (!outArgs.supports(OUT_ARG_DfDp, j).none())
      me_outargs.set_DfDp(j, outArgs.get_DfDp(j));
  for (int i=0; i<num_g; i++) {
    me_outargs.set_g(i, outArgs.get_g(i));
    if (!outArgs.supports(OUT_ARG_DgDx, i).none())
	me_outargs.set_DgDx(i, outArgs.get_DgDx(i));
    if (!outArgs.supports(OUT_ARG_DgDx_dot, i).none())
	me_outargs.set_DgDx(i, outArgs.get_DgDx_dot(i));
    for (int j=0; j<num_p; j++)
      if (!outArgs.supports(OUT_ARG_DgDp, i, j).none())
	me_outargs.set_DgDp(i, j, outArgs.get_DgDp(i,j));
  }

  bool do_quad = false;
  InArgs::sg_const_vector_t x_sg;
  InArgs::sg_const_vector_t x_dot_sg;
  Teuchos::Array<InArgs::sg_const_vector_t> p_sg(num_p);
  OutArgs::sg_vector_t f_sg;
  OutArgs::sg_operator_t W_sg;
  Teuchos::Array<SGDerivative> dfdp_sg(num_p);
  Teuchos::Array<OutArgs::sg_vector_t> g_sg(num_g);
  Teuchos::Array<SGDerivative> dgdx_sg(num_g);
  Teuchos::Array<SGDerivative> dgdx_dot_sg(num_g);
  Teuchos::Array< Teuchos::Array<SGDerivative> > dgdp_sg(num_g);
  TEUCHOS_TEST_FOR_EXCEPTION(inArgs.get_sg_basis() == Teuchos::null, 
		     std::logic_error,
		     "Error!  Stokhos::SGQuadModelEvaluator::evalModel():  " <<
		     "SG basis inArg cannot be null!");
  TEUCHOS_TEST_FOR_EXCEPTION(inArgs.get_sg_quadrature() == Teuchos::null, 
		     std::logic_error,
		     "Error!  Stokhos::SGQuadModelEvaluator::evalModel():  " <<
		     "SG quadrature inArg cannot be null!");
  Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis = 
    inArgs.get_sg_basis();
  Teuchos::RCP< const Stokhos::Quadrature<int,double> > quad = 
    inArgs.get_sg_quadrature();
  if (inArgs.supports(IN_ARG_x_sg)) {
    x_sg = inArgs.get_x_sg();
    if (x_sg != Teuchos::null) {
      do_quad = true;
    }
  }
  if (inArgs.supports(IN_ARG_x_dot_sg)) {
    x_dot_sg = inArgs.get_x_dot_sg();
    if (x_dot_sg != Teuchos::null) {
      do_quad = true;
    }
  }
  for (int i=0; i<num_p; i++) {
    p_sg[i] = inArgs.get_p_sg(i);
    if (p_sg[i] != Teuchos::null) {
      do_quad = true;
    }
  }
  if (outArgs.supports(OUT_ARG_f_sg)) {
    f_sg = outArgs.get_f_sg();
    if (f_sg != Teuchos::null)
      f_sg->init(0.0);
  }
  if (outArgs.supports(OUT_ARG_W_sg)) {
    W_sg = outArgs.get_W_sg();
    if (W_sg != Teuchos::null)
      W_sg->init(0.0);
  }
  for (int i=0; i<num_p; i++) {
    if (!outArgs.supports(OUT_ARG_DfDp_sg, i).none()) {
      dfdp_sg[i] = outArgs.get_DfDp_sg(i);
      if (dfdp_sg[i].getMultiVector() != Teuchos::null)
	dfdp_sg[i].getMultiVector()->init(0.0);
      else if (dfdp_sg[i].getLinearOp() != Teuchos::null)
	dfdp_sg[i].getLinearOp()->init(0.0);
    }
  }
      
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      g_outT = Petra::EpetraVector_To_TpetraVectorNonConst(*g_out, commT);
      app->evaluateResponseT(i, curr_time, x_dotT.get(), x_dotdotT.get(), *xT, sacado_param_vec,
			    *g_outT);
      //convert g_outT to Epetra_Vector g_out
      Petra::TpetraVector_To_EpetraVector(g_outT, *g_out, comm);
    }
  }

  //
  // Stochastic Galerkin
  //
#ifdef ALBANY_SG
  InArgs::sg_const_vector_t x_sg = inArgs.get_x_sg();
  if (x_sg != Teuchos::null) {
    app->init_sg(inArgs.get_sg_basis(),
                 inArgs.get_sg_quadrature(),
                 inArgs.get_sg_expansion(),
                 x_sg->productComm());
    InArgs::sg_const_vector_t x_dot_sg  = Teuchos::null;
    InArgs::sg_const_vector_t x_dot_sg  = Teuchos::null;
    if(num_time_deriv > 0)
      x_dotdot_sg  = inArgs.get_x_dotdot_sg();
    if(num_time_deriv > 1)
      x_dotdot_sg  = inArgs.get_x_dotdot_sg();
    if (x_dot_sg != Teuchos::null || x_dotdot_sg != Teuchos::null) {
      alpha = inArgs.get_alpha();
      beta = inArgs.get_beta();
      curr_time  = inArgs.get_t();
    }
    if (x_dotdot_sg != Teuchos::null) {
      omega = inArgs.get_omega();
    }

    InArgs::sg_const_vector_t epetra_p_sg = inArgs.get_p_sg(0);
    Teuchos::Array<int> p_sg_index;
    for (int i=0; i<num_param_vecs; i++) {
      InArgs::sg_const_vector_t p_sg = inArgs.get_p_sg(i);
      if (p_sg != Teuchos::null) {
        p_sg_index.push_back(i);
        for (int j=0; j<p_sg_vals[i].size(); j++) {
          int num_sg_blocks = p_sg->size();
          p_sg_vals[i][j].reset(app->getStochasticExpansion(), num_sg_blocks);
          p_sg_vals[i][j].copyForWrite();
          for (int l=0; l<num_sg_blocks; l++) {
            p_sg_vals[i][j].fastAccessCoeff(l) = (*p_sg)[l][j];
          }
        }
      }
    }

    OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
    OutArgs::sg_operator_t W_sg = outArgs.get_W_sg();
    bool f_sg_computed = false;

    // W_sg
    if (W_sg != Teuchos::null) {
      Stokhos::VectorOrthogPoly<Epetra_CrsMatrix> W_sg_crs(W_sg->basis(),
                                                           W_sg->map());
      for (int i=0; i<W_sg->size(); i++)
        W_sg_crs.setCoeffPtr(
          i,
          Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_sg->getCoeffPtr(i)));
      app->computeGlobalSGJacobian(alpha, beta, omega, curr_time,
                                   x_dot_sg.get(),  x_dotdot_sg.get(), *x_sg,
                                   sacado_param_vec, p_sg_index, p_sg_vals,
                                   f_sg.get(), W_sg_crs);

void 
MockModelEval_D::
evalModel(const InArgs& inArgs, const OutArgs& outArgs) const
{
  int proc = comm->MyPID();

  // 
  // Deterministic calculation
  //

  // Parse InArgs
  RCP<const Epetra_Vector> p1_in = inArgs.get_p(0);
  if (p1_in == Teuchos::null)
    p1_in = p1_init;
  RCP<const Epetra_Vector> p2_in = inArgs.get_p(1);
  if (p2_in == Teuchos::null)
    p2_in = p2_init;

  RCP<const Epetra_Vector> x_in = inArgs.get_x();

  // Parse OutArgs
  RCP<Epetra_Vector> f_out = outArgs.get_f(); 
  if (f_out != Teuchos::null) {
    double p = (*p1_in)[0];
    double xi = (*p2_in)[0];
    if (proc == 0) {
      double x = (*x_in)[0];
      (*f_out)[0] = x - p + xi;
    }
  }

  RCP<Epetra_CrsMatrix> W_out = 
    Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(outArgs.get_W()); 
  if (W_out != Teuchos::null) {
    if (proc == 0) {
      double val = 1.0;
      int i = 0;
      W_out->ReplaceMyValues(i, 1, &val, &i);
    }
  }

  RCP<Epetra_MultiVector> dfdp1 = outArgs.get_DfDp(0).getMultiVector();
  if (dfdp1 != Teuchos::null) {
    if (proc == 0)
      (*dfdp1)[0][0] = -1.0;
  }
  RCP<Epetra_MultiVector> dfdp2 = outArgs.get_DfDp(1).getMultiVector();
  if (dfdp2 != Teuchos::null) {
    if (proc == 0)
      (*dfdp2)[0][0] = 1.0;
  }

  RCP<Epetra_Vector> g_out = outArgs.get_g(0); 
  if (g_out != Teuchos::null) {
    if (proc == 0) {
      double x = (*x_in)[0];
      (*g_out)[0] = 1.0 / x;
    }
  }
    

  RCP<Epetra_MultiVector> dgdx = outArgs.get_DgDx(0).getMultiVector();
  if (dgdx != Teuchos::null) {
    if (proc == 0) {
      double x = (*x_in)[0];
      (*dgdx)[0][0] = -1.0 / (x*x);
    }
  }

  RCP<Epetra_MultiVector> dgdp1 = outArgs.get_DgDp(0,0).getMultiVector();
  if (dgdp1 != Teuchos::null) {
    if (proc == 0) {
      (*dgdp1)[0][0] = 0.0;
    }
  }
  RCP<Epetra_MultiVector> dgdp2 = outArgs.get_DgDp(0,1).getMultiVector();
  if (dgdp2 != Teuchos::null) {
    if (proc == 0) {
      (*dgdp2)[0][0] = 0.0;
    }
  }

  // 
  // Stochastic calculation
  //

#ifdef Piro_ENABLE_Stokhos
  // Parse InArgs
  RCP<const Stokhos::OrthogPolyBasis<int,double> > basis = 
    inArgs.get_sg_basis();
  RCP<Stokhos::OrthogPolyExpansion<int,double> > expn = 
    inArgs.get_sg_expansion();
  InArgs::sg_const_vector_t x_sg = inArgs.get_x_sg();
  InArgs::sg_const_vector_t p1_sg = inArgs.get_p_sg(0);
  InArgs::sg_const_vector_t p2_sg = inArgs.get_p_sg(1);

  // Parse OutArgs
  OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
  if (f_sg != Teuchos::null && proc == 0) {
    for (int block=0; block<f_sg->size(); block++) {
//.........这里部分代码省略.........