Commit ccd7eecd authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Rosenbrock method extended/modified, Sign in RobinBC changed

parent a81287a5
...@@ -116,6 +116,7 @@ namespace AMDiS { ...@@ -116,6 +116,7 @@ namespace AMDiS {
meshesForProblems[i].insert(meshByName[meshName]); meshesForProblems[i].insert(meshByName[meshName]);
problems[i]->setComponentMesh(j, meshByName[meshName]); problems[i]->setComponentMesh(j, meshByName[meshName]);
MSG((problems[i]->getName() + "->setComponentMesh(" + boost::lexical_cast<std::string>(j) + ") = " + meshName + "\n").c_str());
// feSpace // feSpace
int degree = 1; int degree = 1;
......
...@@ -146,7 +146,12 @@ namespace AMDiS { ...@@ -146,7 +146,12 @@ namespace AMDiS {
addCreator("rowda3", new Rowda3::Creator); addCreator("rowda3", new Rowda3::Creator);
addCreator("ros3p", new Ros3p::Creator); addCreator("ros3p", new Ros3p::Creator);
addCreator("rodasp", new Rodasp::Creator); addCreator("rodasp", new Rodasp::Creator);
addCreator("rosi2pw", new ROSI2PW::Creator); addCreator("rosI2P1", new ROSI2P1::Creator);
addCreator("rosI2P2", new ROSI2P2::Creator);
addCreator("rosI2Pw", new ROSI2Pw::Creator);
addCreator("rosI2PW", new ROSI2PW::Creator);
addCreator("ros3Pw", new Ros3Pw::Creator);
addCreator("ros34PW2", new Ros34PW2::Creator);
} }
......
...@@ -209,6 +209,16 @@ struct Vec3WorldVec : public TertiaryAbstractFunction<WorldVector<T>,T,T,T> ...@@ -209,6 +209,16 @@ struct Vec3WorldVec : public TertiaryAbstractFunction<WorldVector<T>,T,T,T>
} }
}; };
struct Component : public AbstractFunction<double, WorldVector<double> >
{
Component(int comp_) : comp(comp_) {}
double operator()(const WorldVector<double> &x) const {
return x[comp];
}
private:
int comp;
};
struct FadeOut : public TertiaryAbstractFunction<double, double, double ,double> struct FadeOut : public TertiaryAbstractFunction<double, double, double ,double>
{ {
double operator()(const double &v, const double &dist, const double &mean) const { double operator()(const double &v, const double &dist, const double &mean) const {
......
...@@ -665,7 +665,7 @@ namespace AMDiS { ...@@ -665,7 +665,7 @@ namespace AMDiS {
return new ElInfo3d(this); return new ElInfo3d(this);
break; break;
default: default:
ERROR_EXIT("invalid dim\n"); ERROR_EXIT("invalid dim [%d]\n",dim);
return NULL; return NULL;
} }
} }
......
...@@ -68,6 +68,8 @@ namespace AMDiS { ...@@ -68,6 +68,8 @@ namespace AMDiS {
(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]); (*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);
delete alphaOp; delete alphaOp;
WARNING("Sign of alpha changed in RobinBC!\n");
} }
} }
...@@ -119,6 +121,7 @@ namespace AMDiS { ...@@ -119,6 +121,7 @@ namespace AMDiS {
(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]); (*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);
delete alphaOp; delete alphaOp;
WARNING("Sign of alpha changed in RobinBC!\n");
} }
} }
...@@ -159,6 +162,7 @@ namespace AMDiS { ...@@ -159,6 +162,7 @@ namespace AMDiS {
robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT); robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
for (int i = 0; i < dim + 1; i++) for (int i = 0; i < dim + 1; i++)
(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]); (*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);
WARNING("Sign of alpha changed in RobinBC!\n");
} }
} }
...@@ -193,7 +197,7 @@ namespace AMDiS { ...@@ -193,7 +197,7 @@ namespace AMDiS {
for (int i = 0; i < dim + 1; i++) for (int i = 0; i < dim + 1; i++)
if (elInfo->getBoundary(i) == boundaryType) if (elInfo->getBoundary(i) == boundaryType)
matrix->assemble(-1.0, elInfo, localBound, (*robinOperators)[i]); matrix->assemble(1.0, elInfo, localBound, (*robinOperators)[i]);
} }
} }
...@@ -262,7 +266,7 @@ namespace AMDiS { ...@@ -262,7 +266,7 @@ namespace AMDiS {
f = 0.0; f = 0.0;
if (robinOperators) if (robinOperators)
(*robinOperators)[face]->evalZeroOrder(nPoints, uhAtQp, grdUhAtQp, D2UhAtQp, f, 1.0); (*robinOperators)[face]->evalZeroOrder(nPoints, uhAtQp, grdUhAtQp, D2UhAtQp, f, -1.0);
std::vector<WorldVector<double> > grdUh(nPoints); std::vector<WorldVector<double> > grdUh(nPoints);
std::vector<WorldVector<double> > A_grdUh(nPoints); std::vector<WorldVector<double> > A_grdUh(nPoints);
......
...@@ -34,7 +34,7 @@ namespace AMDiS { ...@@ -34,7 +34,7 @@ namespace AMDiS {
* *
* \brief * \brief
* Sub class of BoundaryCondition. Implements Robin and Neumann boundary conditions. * Sub class of BoundaryCondition. Implements Robin and Neumann boundary conditions.
* The flux in normal direction is given by \f$ j = j_0 + \alpha u\f$ where * The flux in normal direction is given by \f$ A\nabla u\cdot \nu\phi =: j = j_0 - \alpha u \f$ where
* \f$ j_0 \f$ and \f$ alpha \f$ are functions evaluated at world coordinates * \f$ j_0 \f$ and \f$ alpha \f$ are functions evaluated at world coordinates
* and \f$ u \f$ is the problem solution. * and \f$ u \f$ is the problem solution.
*/ */
......
...@@ -45,7 +45,7 @@ void RobinBC<T>::fillLocalBC(DOFMatrix* matrix, ...@@ -45,7 +45,7 @@ void RobinBC<T>::fillLocalBC(DOFMatrix* matrix,
if(robinOperators) { if(robinOperators) {
for(i=0; i < dim+1; i++) { for(i=0; i < dim+1; i++) {
if(elInfo->getBoundOfBoundary(i) == boundaryType) { if(elInfo->getBoundOfBoundary(i) == boundaryType) {
matrix->assemble(-1.0, elInfo, localBound, (*robinOperators)[i]); matrix->assemble(1.0, elInfo, localBound, (*robinOperators)[i]);
} }
} }
} }
......
...@@ -46,17 +46,30 @@ namespace AMDiS { ...@@ -46,17 +46,30 @@ namespace AMDiS {
gamma = 1.707106781186547; gamma = 1.707106781186547;
createData(); createData();
// b1(2) = 0.5
// b1(1) = 1-b1(2) = 0.5
// alpha(2,1) = 1/(2*b1(2)) = 1
// gamma = 1+1/sqrt(2) = 1.707...
// gamma(2,1) = -gamma/b1(2) = -3.41421356237309
// b2(1) = 1
// b2(2) = 0
a[0][0] = 0.0;
a[1][0] = 5.857864376269050e-01; a[1][0] = 5.857864376269050e-01;
a[1][1] = a[1][0];
c[0][0] = 5.857864376269050e-01; c[0][0] = gamma;
c[1][0] = 1.171572875253810e+00; c[1][0] = -1.171572875253810e+00;
c[1][1] = 5.857864376269050e-01; c[1][1] = c[1][0] + gamma;
m1[0] = 8.786796564403575e-01; m1[0] = 8.786796564403575e-01;
m1[1] = 2.928932188134525e-01; m1[1] = 2.928932188134525e-01;
m2[0] = 5.857864376269050e-01; m2[0] = 5.857864376269050e-01;
m2[1] = 0.0;
MSG("Rosenbrock scheme Ros2\n");
} }
...@@ -68,15 +81,18 @@ namespace AMDiS { ...@@ -68,15 +81,18 @@ namespace AMDiS {
createData(); createData();
a[0][0] = 0.0;
a[1][0] = 1.605996252195329e+00; a[1][0] = 1.605996252195329e+00;
a[1][1] = 0.7;
a[2][0] = 1.605996252195329e+00; a[2][0] = 1.605996252195329e+00;
a[2][2] = 0.7;
c[0][0] = 2.294280360279042e+00; c[0][0] = 0.435866521508459;
c[1][0] = -8.874044410657833e-01; c[1][0] = 8.874044410657833e-01;
c[1][1] = 2.294280360279042e+00; c[1][1] = 0.604455284065559;
c[2][0] = -2.398747971635036e+01; c[2][0] = 2.398747971635036e+01;
c[2][1] = -5.263722371562129e+00; c[2][1] = 5.263722371562129e+00;
c[2][2] = 2.294280360279042e+00; c[2][2] = 6.37978879934488;
m1[0] = 2.236727045296590e+00; m1[0] = 2.236727045296590e+00;
m1[1] = 2.250067730969644e+00; m1[1] = 2.250067730969644e+00;
...@@ -84,6 +100,8 @@ namespace AMDiS { ...@@ -84,6 +100,8 @@ namespace AMDiS {
m2[0] = 2.059356167645940e+00; m2[0] = 2.059356167645940e+00;
m2[1] = 1.694014319346528e-01; m2[1] = 1.694014319346528e-01;
MSG("Rosenbrock scheme Rowda3\n");
} }
...@@ -95,23 +113,29 @@ namespace AMDiS { ...@@ -95,23 +113,29 @@ namespace AMDiS {
createData(); createData();
a[1][0] = 1.267949192431123e+00; a[0][0] = 0.0;
a[2][0] = 1.267949192431123e+00; a[1][0] = 1.26794919243112;
a[1][1] = 1.0;
a[2][0] = 1.26794919243112;
a[2][1] = 0.0;
a[2][2] = 1.0;
c[0][0] = 1.267949192431123e+00; c[0][0] = 0.788675134594813;
c[1][0] = 1.607695154586739e+00; c[1][0] = -1.60769515458674;
c[1][1] = 1.267949192431123e+00; c[1][1] = -0.211324865405187;
c[2][0] = 3.462101615137755e+00; c[2][0] = -3.46410161513775;
c[2][1] = 1.732050807568877e+00; c[2][1] = -1.73205080756888;
c[2][2] = 1.267949192431123e+00; c[2][2] = -1.07735026918963;
m1[0] = 2.0; m1[0] = 2.0;
m1[1] = 5.773502691896258e-01; m1[1] = 0.577350269189625;
m1[2] = 4.226497308103742e-01; m1[2] = 0.422649730810374;
m2[0] = 2.113248654051871e+00; m2[0] = 2.11324865405187;
m2[1] = 1.0; m2[1] = 1.0;
m2[2] = 4.226497308103742e-01; m2[2] = 0.422649730810374;
MSG("Rosenbrock scheme Ros3p\n");
} }
...@@ -123,43 +147,51 @@ namespace AMDiS { ...@@ -123,43 +147,51 @@ namespace AMDiS {
createData(); createData();
a[0][0] = 0.0;
a[1][0] = 3.0; a[1][0] = 3.0;
// a[1][1] =
a[2][0] = 1.831036793486759e+00; a[2][0] = 1.831036793486759e+00;
a[2][1] = 4.955183967433795e-01; a[2][1] = 4.955183967433795e-01;
// a[2][2] =
a[3][0] = 2.304376582692669e+00; a[3][0] = 2.304376582692669e+00;
a[3][1] = -5.249275245743001e-02; a[3][1] = -5.249275245743001e-02;
a[3][2] = -1.176798761832782e+00; a[3][2] = -1.176798761832782e+00;
// a[3][3] =
a[4][0] = -7.170454962423025e+00; a[4][0] = -7.170454962423025e+00;
a[4][1] = -4.741636671481786e+00; a[4][1] = -4.741636671481786e+00;
a[4][2] = -1.631002631330971e+01; a[4][2] = -1.631002631330971e+01;
a[4][3] = -1.062004044111401e+00; a[4][3] = -1.062004044111401e+00;
// a[4][4] =
a[5][0] = -7.170454962423025e+00; a[5][0] = -7.170454962423025e+00;
a[5][1] = -4.741636671481785e+00; a[5][1] = -4.741636671481785e+00;
a[5][2] = -1.631002631330971e+01; a[5][2] = -1.631002631330971e+01;
a[5][3] = -1.062004044111401e+00; a[5][3] = -1.062004044111401e+00;
a[5][4] = 1.0; a[5][4] = 1.0;
a[5][5] = 1.0;
c[0][0] = 4.0; c[0][0] = gamma;
c[1][0] = 1.200000000000000e+01;
c[1][1] = 4.0; c[1][0] = -1.200000000000000e+01;
c[2][0] = 8.791795173947035e+00; // c[1][1] = 4.0;
c[2][1] = 2.207865586973518e+00; c[2][0] = -8.791795173947035e+00;
c[2][2] = 4.0; c[2][1] = -2.207865586973518e+00;
c[3][0] = -1.081793056857153e+01; // c[2][2] = 4.0;
c[3][1] = -6.780270611428266e+00; c[3][0] = 1.081793056857153e+01;
c[3][2] = -1.953485944642410e+01; c[3][1] = 6.780270611428266e+00;
c[3][3] = 4.0; c[3][2] = 1.953485944642410e+01;
c[4][0] = -3.419095006749677e+01; // c[3][3] = 4.0;
c[4][1] = -1.549671153725963e+01; c[4][0] = 3.419095006749677e+01;
c[4][2] = -5.474760875964130e+01; c[4][1] = 1.549671153725963e+01;
c[4][3] = -1.416005392148534e+01; c[4][2] = 5.474760875964130e+01;
c[4][4] = 4.0; c[4][3] = 1.416005392148534e+01;
c[5][0] = -3.462605830930533e+01; // c[4][4] = 4.0;
c[5][1] = -1.530084976114473e+01; c[5][0] = 3.462605830930533e+01;
c[5][2] = -5.699955578662667e+01; c[5][1] = 1.530084976114473e+01;
c[5][3] = -1.840807009793095e+01; c[5][2] = 5.699955578662667e+01;
c[5][4] = 5.714285714285717e+00; c[5][3] = 1.840807009793095e+01;
c[5][5] = 4.0; c[5][4] = -5.714285714285717e+00;
// c[5][5] = 4.0;
m1[0] = -7.170454962423026e+00; m1[0] = -7.170454962423026e+00;
m1[1] = -4.741636671481786e+00; m1[1] = -4.741636671481786e+00;
...@@ -173,8 +205,142 @@ namespace AMDiS { ...@@ -173,8 +205,142 @@ namespace AMDiS {
m2[2] = -1.631002631330971e+01; m2[2] = -1.631002631330971e+01;
m2[3] = -1.062004044111401e+00; m2[3] = -1.062004044111401e+00;
m2[4] = 1.0; m2[4] = 1.0;
WARNING("You have to fill the diagonals of the coeffcients to use time derivatives of your operators!\n");
MSG("Rosenbrock scheme Rodasp\n");
}
ROSI2P1::ROSI2P1()
{
order = 3;
stages = 4;
gamma = 4.35866521508459e-01;
createData();
a[0][0] = 0.0;
a[1][0] = 1.14714018013952;
a[1][1] = 0.5;
a[2][0] = 1.32930341703716;
a[2][1] = 0.442124760965983;
a[2][2] = 0.75;
a[3][0] = -2.19977665049954;
a[3][1] = 4.55821087656518;
a[3][2] = -1.37361554490645;
a[3][3] = 1.0;
c[0][0] = gamma;
c[1][0] = -0.263186118578107;
c[1][1] = 0.385866521508459;
c[2][0] = -3.35635061779949;
c[2][1] = 0.334203214637756;
c[2][2] = -0.145563307177156;
c[3][0] = -2.6183067448488;
c[3][1] = -1.08994123815716;
c[3][2] = -1.71836543021444;
c[3][3] = -0.135847884055848;
m1[0] = -1.24823690752574;
m1[1] = 3.9534178869996;
m1[2] = -1.21522771421847;
m1[3] = 1.16541744745931;
m2[0] = 0.919998349829964;
m2[1] = 1.77038079363523;
m2[2] = 0.257316038155499;
m2[3] = 0.343556220548095;
MSG("Rosenbrock scheme ROSI2P1\n");
} }
ROSI2P2::ROSI2P2()
{
order = 3;
stages = 4;
gamma = 4.3586652150845900e-01;
createData();
a[0][0] = 0.0;
a[1][0] = 1.14714018013952;
a[1][1] = 0.5;
a[2][0] = -0.79265181178863;
a[2][1] = 3.48693217206767;
a[2][2] = 1.0;
a[3][0] = -0.79265181178863;
a[3][1] = 3.48693217206767;
a[3][2] = 0.0;
a[3][3] = 1.0;
c[0][0] = gamma;
c[1][0] = -0.263186118578107;
c[1][1] = 0.385866521508459;
c[2][0] = -1.40507848232679;
c[2][1] = 6.18104102134041;
c[2][2] = 1.20849664917601;
c[3][0] = 4.99718323859189;
c[3][1] = -6.54597265243973;
c[3][2] = -0.539706236424999;
c[3][3] = 0.0;
m1[0] = -0.792651811788631;
m1[1] = 3.48693217206767;
m1[2] = 0.0;
m1[3] = 1.0;
m2[0] = 6.55610198824936;
m2[1] = -5.94329934171132;
m2[2] = 0.360559439940373;
m2[3] = -3.34373891242489;
MSG("Rosenbrock scheme ROSI2P2\n");
}
ROSI2Pw::ROSI2Pw()
{
order = 3;
stages = 4;
gamma = 4.3586652150845e-01;
createData();
a[0][0] = 0.0;
a[1][0] = 2.0;
a[1][1] = 0.871733043016918;
a[2][0] = 1.63034046718537;
a[2][1] = -0.0903698030239437;
a[2][2] = 0.75;
a[3][0] = 3.28685750853379;
a[3][1] = 15.4503745896552;
a[3][2] = -15.0445558288081;
a[3][3] = 1.0;
c[0][0] = gamma;
c[1][0] = -4.58856072055827;
c[1][1] = -0.435866521508468;
c[2][0] = -4.56739138878308;
c[2][1] = -0.0683107605436897;
c[2][2] = -0.418867127163069;
c[3][0] = -2.99296365068231;
c[3][1] = -42.6472578877056;
c[3][2] = 43.6510246478391;
c[3][3] = 0.0;
m1[0] = 3.28685750853379;
m1[1] = 15.4503745896552;
m1[2] = -15.0445558288081;
m1[3] = 1.0;
m2[0] = 3.3904377475357;
m2[1] = 5.865078412615;
m2[2] = -4.96246395067988;
m2[3] = 0.260825116305751;
MSG("Rosenbrock scheme ROSI2Pw\n");
}
ROSI2PW::ROSI2PW() ROSI2PW::ROSI2PW()
{ {
...@@ -183,29 +349,124 @@ namespace AMDiS { ...@@ -183,29 +349,124 @@ namespace AMDiS {
gamma = 4.3586652150845e-01; gamma = 4.3586652150845e-01;
createData(); createData();
a[0][0] = 0.0;
a[1][0] = 2.0;
a[1][1] = 0.871733043016918;
a[2][0] = -5.501959790112310;
a[2][1] = -1.833986596704078;
a[2][2] = -1.59874671679705;
a[3][0] = 4.338560720558247;
a[3][1] = 1.147140180139553;
a[3][2] = -0.059559354637499;
a[3][3] = 1.0;
c[0][0] = gamma;
c[1][0] = -4.58856072055827;
c[1][1] = -0.435866521508468;
c[2][0] = 48.3965596116246;
c[2][1] = 16.132186537208;
c[2][2] = 6.56544000523295;
c[3][0] = -2.31911621201727;
c[3][1] = -1.14714018013959;
c[3][2] = -0.0745075552140216;
c[3][3] = 0.0;
m1[0] = 4.33856072055832;
m1[1] = 1.14714018013958;
m1[2] = -0.0595593546375007;
m1[3] = 1.0;
m2[0] = 3.31383144448529;
m2[1] = 1.14714018013955;
m2[2] = 0.00847038665840751;
m2[3] = 0.260825116305751;
MSG("Rosenbrock scheme ROSI2PW\n");
}
Ros3Pw::Ros3Pw()
{
// index1 | index2 | pdes | R(infty) | stiffly acc.
// -------+--------+------+----------+-------------
// 1 | 0 | 1 | 0.73 | 0
order = 3;
stages = 3;
gamma = 7.8867513459481287e-01;
createData();
a[0][0] = 0.0;
a[1][0] = 2.0;
a[1][1] = 1.57735026918963;
a[2][0] = 0.633974596215561;
a[2][1] = 0.0;
a[2][2] = 0.5;
c[0][0] = 0.788675134594813;
c[1][0] = -2.53589838486225;
c[1][1] = -0.788675134594813;
c[2][0] = -1.62740473580836;