diff --git a/test/decoupled/1p/test_1p-reference.vtu b/test/decoupled/1p/test_1p-reference.vtu
index aafd1798cd40d6f2dc8d7762ae7653119f4d5e37..40e13bb5341cabfc62edf2d95d04672eec5586e9 100644
--- a/test/decoupled/1p/test_1p-reference.vtu
+++ b/test/decoupled/1p/test_1p-reference.vtu
@@ -4,30 +4,30 @@
<Piece NumberOfCells="64" NumberOfPoints="81">
<CellData Scalars="pressure" Vectors="velocity">
<DataArray type="Float32" Name="pressure" NumberOfComponents="1" format="ascii">
- 0.193321 0.276067 0.279825 0.244879 0.188075 0.119172 0.0457015 -0.0186704 0.276067 0.644802 0.738844 0.676876
- 0.535036 0.350018 0.1468 -0.0333657 0.279825 0.738844 0.958362 0.951347 0.795798 0.553011 0.267066 0.0014723
- 0.244879 0.676876 0.951347 1.02608 0.925178 0.697408 0.391748 0.0757681 0.188075 0.535036 0.795798 0.925178
- 0.906034 0.749354 0.483548 0.155017 0.119172 0.350018 0.553011 0.697408 0.749354 0.685392 0.500228 0.203568
- 0.0457015 0.1468 0.267066 0.391748 0.483548 0.500228 0.411616 0.196857 -0.0186704 -0.0333657 0.0014723 0.0757681
- 0.155017 0.203568 0.196857 0.113592
+ 0.134323 0.183052 0.203147 0.19449 0.161069 0.109222 0.0462204 -0.0147278 0.183052 0.395195 0.512464 0.526295
+ 0.45381 0.319656 0.147713 -0.0222035 0.203147 0.512464 0.702329 0.758698 0.686768 0.512331 0.269591 0.0171451
+ 0.19449 0.526295 0.758698 0.86232 0.826021 0.661429 0.398211 0.0935879 0.161069 0.45381 0.686768 0.826021
+ 0.844359 0.731349 0.496948 0.174174 0.109222 0.319656 0.512331 0.661429 0.731349 0.690355 0.521748 0.223249
+ 0.0462204 0.147713 0.269591 0.398211 0.496948 0.521748 0.437118 0.21426 -0.0147278 -0.0222035 0.0171451 0.0935879
+ 0.174174 0.223249 0.21426 0.122997
</DataArray>
<DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="ascii">
- -0.939713 -0.939713 0 -0.0349131 -3.31549 0 0.00491799 -3.91811 0 0.00769967 -3.61336 0
- 0.00662874 -2.8572 0 0.0052328 -1.8614 0 0.00379149 -0.765481 0 0.000586072 0.207224 0
- -3.31549 -0.0349131 0 -0.926479 -0.926479 0 -0.0340556 -1.99642 0 0.12719 -2.38781 0
- 0.13192 -2.18805 0 0.108957 -1.61534 0 0.0789943 -0.840733 0 0.0178877 -0.0774749 0
- -3.91811 0.00491799 0 -1.99642 -0.0340556 0 -0.425431 -0.425431 0 0.220112 -0.925378 0
- 0.377005 -1.19288 0 0.363898 -1.15186 0 0.28622 -0.853655 0 0.108382 -0.392925 0
- -3.61336 0.00769967 0 -2.38781 0.12719 0 -0.925378 0.220112 0 0.0523886 0.0523886 0
- 0.496349 -0.27491 0 0.616533 -0.559225 0 0.560833 -0.671487 0 0.335963 -0.504454 0
- -2.8572 0.00662874 0 -2.18805 0.13192 0 -1.19288 0.377005 0 -0.27491 0.496349 0
- 0.352001 0.352001 0 0.678664 0.0288108 0 0.771867 -0.293473 0 0.678006 -0.376017 0
- -1.8614 0.0052328 0 -1.61534 0.108957 0 -1.15186 0.363898 0 -0.559225 0.616533 0
- 0.0288108 0.678664 0 0.498749 0.498749 0 0.80527 0.167795 0 0.986332 -0.108892 0
- -0.765481 0.00379149 0 -0.840733 0.0789943 0 -0.853655 0.28622 0 -0.671487 0.560833 0
- -0.293473 0.771867 0 0.167795 0.80527 0 0.607349 0.607349 0 1.04537 0.205684 0
- 0.207224 0.000586072 0 -0.0774749 0.0178877 0 -0.392925 0.108382 0 -0.504454 0.335963 0
- -0.376017 0.678006 0 -0.108892 0.986332 0 0.205684 1.04537 0 0.621519 0.621519 0
+ -0.635385 -0.635385 0 -0.105666 -2.08192 0 -0.00674492 -2.75246 0 0.00416985 -2.82553 0
+ 0.00553863 -2.42955 0 0.00504952 -1.70147 0 0.00399201 -0.771175 0 0.00103779 0.147072 0
+ -2.08192 -0.105666 0 -0.99857 -0.99857 0 -0.249458 -1.74836 0 0.0305497 -2.08621 0
+ 0.0993801 -1.99544 0 0.1024 -1.55224 0 0.0832657 -0.866472 0 0.0280078 -0.121576 0
+ -2.75246 -0.00674492 0 -1.74836 -0.249458 0 -0.67148 -0.67148 0 -0.0154582 -1.0601 0
+ 0.26249 -1.26451 0 0.332161 -1.21177 0 0.298072 -0.911135 0 0.144846 -0.426571 0
+ -2.82553 0.00416985 0 -2.08621 0.0305497 0 -1.0601 -0.0154582 0 -0.201927 -0.201927 0
+ 0.321234 -0.471503 0 0.552789 -0.693519 0 0.57905 -0.754513 0 0.409333 -0.53998 0
+ -2.42955 0.00553863 0 -1.99544 0.0993801 0 -1.26451 0.26249 0 -0.471503 0.321234 0
+ 0.180515 0.180515 0 0.602776 -0.100857 0 0.79572 -0.373355 0 0.789052 -0.409145 0
+ -1.70147 0.00504952 0 -1.55224 0.1024 0 -1.21177 0.332161 0 -0.693519 0.552789 0
+ -0.100857 0.602776 0 0.435841 0.435841 0 0.836432 0.130326 0 1.12462 -0.120997 0
+ -0.771175 0.00399201 0 -0.866472 0.0832657 0 -0.911135 0.298072 0 -0.754513 0.57905 0
+ -0.373355 0.79572 0 0.130326 0.836432 0 0.643578 0.643578 0 1.18237 0.232018 0
+ 0.147072 0.00103779 0 -0.121576 0.0280078 0 -0.426571 0.144846 0 -0.53998 0.409333 0
+ -0.409145 0.789052 0 -0.120997 1.12462 0 0.232018 1.18237 0 0.701106 0.701106 0
</DataArray>
</CellData>
<Points>
diff --git a/test/decoupled/1p/test_1p_problem.hh b/test/decoupled/1p/test_1p_problem.hh
index 5039bd570bfde17024ee0f5842bef5534dc73ded..1049ce9b501fffce0ef93ed8c794ea91155f3bef 100644
--- a/test/decoupled/1p/test_1p_problem.hh
+++ b/test/decoupled/1p/test_1p_problem.hh
@@ -116,6 +116,7 @@ class TestProblemOneP: public DiffusionProblem1P<TypeTag >
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::Intersection Intersection;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
@@ -175,7 +176,7 @@ public:
*
* This problem assumes a temperature of 10 degrees Celsius.
*/
- Scalar temperature(const Element& element) const
+ Scalar temperatureAtPos(const GlobalPosition& globalPos) const
{
return 273.15 + 10; // -> 10°C
}
@@ -183,26 +184,17 @@ public:
// \}
//! Returns the reference pressure for evaluation of constitutive relations
- Scalar referencePressure(const Element& element) const
+ Scalar referencePressureAtPos(const GlobalPosition& globalPos) const
{
return 1e5; // -> 10°C
}
//!source term [kg/(m^3 s)]
- void sourceAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
+ void source(PrimaryVariables &values, const Element& element) const
{
- double pi = 4.0*atan(1.0);
- double rt = globalPos[0]*globalPos[0]+globalPos[1]*globalPos[1];
- double ux = pi*cos(pi*globalPos[0])*sin(pi*globalPos[1]);
- double uy = pi*cos(pi*globalPos[1])*sin(pi*globalPos[0]);
- double kxx = (delta_*globalPos[0]*globalPos[0] + globalPos[1]*globalPos[1])/rt;
- double kxy = -(1.0 - delta_)*globalPos[0]*globalPos[1]/rt;
- double kyy = (globalPos[0]*globalPos[0] + delta_*globalPos[1]*globalPos[1])/rt;
- double f0 = sin(pi*globalPos[0])*sin(pi*globalPos[1])*pi*pi*(1.0 + delta_)*(globalPos[0]*globalPos[0] + globalPos[1]*globalPos[1])
- + cos(pi*globalPos[0])*sin(pi*globalPos[1])*pi*(1.0 - 3.0*delta_)*globalPos[0]
- + cos(pi*globalPos[1])*sin(pi*globalPos[0])*pi*(1.0 - 3.0*delta_)*globalPos[1]
- + cos(pi*globalPos[1])*cos(pi*globalPos[0])*2.0*pi*pi*(1.0 - delta_)*globalPos[0]*globalPos[1];
- values = (f0 + 2.0*(globalPos[0]*(kxx*ux + kxy*uy) + globalPos[1]*(kxy*ux + kyy*uy)))/rt;
+ values = 0;
+
+ values = integratedSource_(element, 4);
}
/*!
@@ -248,6 +240,55 @@ private:
return grad;
}
+ Scalar integratedSource_(const Element& element, int integrationPoints) const
+ {
+ Scalar source = 0.;
+ LocalPosition localPos(0.0);
+ GlobalPosition globalPos(0.0);
+ Scalar halfInterval = 1.0/double(integrationPoints)/2.;
+ for (int i = 1; i <= integrationPoints; i++)
+ {
+ for (int j = 1; j <= integrationPoints; j++)
+ {
+ localPos[0] = double(i)/double(integrationPoints) - halfInterval;
+ localPos[1] = double(j)/double(integrationPoints) - halfInterval;
+ globalPos = element.geometry().global(localPos);
+ source += 1./(integrationPoints*integrationPoints) * evaluateSource_(globalPos);
+ }
+ }
+
+ return source;
+ }
+
+ Scalar evaluateSource_(const GlobalPosition& globalPos) const
+ {
+ Scalar temp = temperatureAtPos(globalPos);
+ Scalar referencePress = referencePressureAtPos(globalPos);
+
+ Scalar pi = 4.0 * atan(1.0);
+ Scalar x = globalPos[0];
+ Scalar y = globalPos[1];
+
+ Scalar dpdx = pi * cos(pi * x) * sin(pi * y);
+ Scalar dpdy = pi * sin(pi * x) * cos(pi * y);
+ Scalar dppdxx = -pi * pi * sin(pi * x) * sin(pi * y);
+ Scalar dppdxy = pi * pi * cos(pi * x) * cos(pi * y);
+ Scalar dppdyx = dppdxy;
+ Scalar dppdyy = dppdxx;
+ Scalar kxx = (delta_* x*x + y*y)/(x*x + y*y);
+ Scalar kxy = -(1.0 - delta_) * x * y / (x*x + y*y);
+ Scalar kyy = (x*x + delta_*y*y)/(x*x + y*y);
+ Scalar dkxxdx = 2 * x * y*y * (delta_ - 1.0)/((x*x + y*y) * (x*x + y*y));
+ Scalar dkyydy = 2 * x*x * y * (delta_ - 1.0)/((x*x + y*y) * (x*x + y*y));
+ Scalar dkxydx = (1.0 - delta_) * y * (x*x - y*y) /((x*x + y*y) * (x*x + y*y));
+ Scalar dkxydy = (1.0 - delta_) * x * (y*y - x*x) /((x*x + y*y) * (x*x + y*y));
+
+ Scalar fx = dkxxdx * dpdx + kxx * dppdxx + dkxydx * dpdy + kxy * dppdyx;
+ Scalar fy = dkxydy * dpdx + kxy * dppdxy + dkyydy * dpdy + kyy * dppdyy;
+
+ return -(fx + fy) / Fluid::viscosity(temp, referencePress) * Fluid::density(temp, referencePress);
+ }
+
double delta_;
Dumux::FVVelocity<TypeTag, typename GET_PROP_TYPE(TypeTag, Velocity) > velocity_;
};