PlayaPoissonBoltzmannOp.cpp
00001
00002
00003
00004
00005 #include "PlayaPoissonBoltzmannOp.hpp"
00006 #include "PlayaTabs.hpp"
00007
00008 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00009 #include "PlayaVectorImpl.hpp"
00010 #include "PlayaLinearOperatorImpl.hpp"
00011 #endif
00012
00013 using namespace Teuchos;
00014
00015 namespace Playa
00016 {
00017
00018 PoissonBoltzmannOp::PoissonBoltzmannOp(int nLocal, const VectorType<double>& vecType)
00019 : NonlinearOperatorBase<double>(), J_(nLocal, vecType), importer_(),
00020 uLeftBC_(0.0), uRightBC_(2.0*log(cosh(1.0/sqrt(2.0))))
00021 {
00022 setDomainAndRange(J_.domain(), J_.range());
00023
00024 int rank = MPIComm::world().getRank();
00025 int nProc = MPIComm::world().getNProc();
00026 if (nProc > 1)
00027 {
00028 Array<int> ghosts;
00029 int low = J_.domain().baseGlobalNaturalIndex();
00030 int high = low + J_.domain().numLocalElements();
00031 if (rank != nProc - 1)
00032 {
00033 ghosts.append(high);
00034 }
00035 if (rank != 0)
00036 {
00037 ghosts.append(low-1);
00038 }
00039
00040 importer_ = vecType.createGhostImporter(J_.domain(), ghosts.size(), &(ghosts[0]));
00041 }
00042 else
00043 {
00044 importer_ = vecType.createGhostImporter(J_.domain(), 0, 0);
00045 }
00046 }
00047
00048 Vector<double> PoissonBoltzmannOp::getInitialGuess() const
00049 {
00050 Vector<double> rtn = J_.domain().createMember();
00051
00052 rtn.setToConstant(0.5);
00053
00054 return rtn;
00055 }
00056
00057
00058 LinearOperator<double>
00059 PoissonBoltzmannOp::computeJacobianAndFunction(Vector<double>& functionValue) const
00060 {
00061 Tabs tab;
00062 Out::root() << tab << "in PBOp::computeJacAndVec" << std::endl;
00063 Vector<double> x = currentEvalPt() ;
00064 Out::root() << tab << "eval pt = " << std::endl;
00065 Out::os() << x << std::endl;
00066 J_.setEvalPoint(x);
00067
00068
00069 RCP<GhostView<double> > u;
00070 Out::root() << tab << "importing view" << std::endl;
00071 importer_->importView(currentEvalPt(), u);
00072 Out::root() << tab << "done importing view" << std::endl;
00073 int low = J_.domain().baseGlobalNaturalIndex();
00074 int high = low + J_.domain().numLocalElements();
00075 Out::os() << tab << "my indices are: " << low << ", " << high-1 << std::endl;
00076
00077 functionValue = J_.range().createMember();
00078 double h= J_.h();
00079
00080 for (int r=low; r<high; r++)
00081 {
00082 Tabs tab1;
00083 double u_i = u->getElement(r);
00084 double f = 0.0;
00085 if (r==0)
00086 {
00087 f = u_i - uLeftBC_;
00088 }
00089 else if (r==J_.domain().dim()-1)
00090 {
00091 f = u_i - uRightBC_;
00092 }
00093 else
00094 {
00095 double u_plus = u->getElement(r+1);
00096 double u_minus = u->getElement(r-1);
00097 f = (u_plus + u_minus - 2.0*u_i)/h/h - exp(-u_i);
00098 }
00099 functionValue[r-low] = f;
00100 }
00101
00102 Out::root() << tab << "done PBOp::computeJacAndVec" << std::endl;
00103 return J_.getOp();
00104 }
00105
00106 Vector<double> PoissonBoltzmannOp::exactSoln() const
00107 {
00108 Tabs tab;
00109 Out::root() << tab << "in PBOp::exactSoln" << std::endl;
00110 Vector<double> rtn = J_.domain().createMember();
00111
00112 int low = J_.domain().baseGlobalNaturalIndex();
00113 int high = low + J_.domain().numLocalElements();
00114
00115 double h= J_.h();
00116
00117 for (int r=low; r<high; r++)
00118 {
00119 double x = r*h;
00120 double u = 2.0*log(cosh(x/sqrt(2.0)));
00121 rtn[r-low] = u;
00122 }
00123
00124 Out::os() << tab << "done PBOp::exactSoln" << std::endl;
00125 return rtn;
00126 }
00127
00128 }