From 6f005bcfa811d7dd163c31353d414bf9d4bccdfa Mon Sep 17 00:00:00 2001
From: Ondrej Szekely <ondra.szekely@gmail.com>
Date: Wed, 29 Jul 2015 16:06:25 +0200
Subject: [PATCH] THE LAST FIXING!!!!
---
.../tnl-run-mean-curvature-flow-eoc-test | 26 +++++++++----------
src/functors/tnlSinWaveFunction_impl.h | 12 ++++++---
.../tnlFiniteVolumeNonlinearOperator_impl.h | 12 ++++-----
.../tnlFiniteVolumeOperatorQ_impl.h | 24 ++++++++---------
4 files changed, 38 insertions(+), 36 deletions(-)
diff --git a/examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test b/examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test
index aa151c5a8a..c4b9057e6c 100755
--- a/examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test
+++ b/examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test
@@ -1,23 +1,22 @@
#!/bin/bash
device="host"
-dimensions="2D"
-sizes2D="16 32"
-testFunctions="sin-bumps"
-snapshotPeriod=0.1
-finalTime=0.3
+dimensions="3D"
+sizes3D="16 32 64"
+testFunctions="exp-bump sin-wave sin-bumps"
+snapshotPeriod=0.001
+finalTime=0.01
timeDependence="cosine"
solverName="mean-curvature-flow-eoc"
#solverName="gdb --args tnl-heat-equation-eoc-test-dbg"
-#
setupTestFunction()
{
testFunction=$1
# if test x${testFunction} = "xexp-bump";
# then
- origin=0
- proportions=1
+ origin=2
+ proportions=2
amplitude=1.0
waveLength=1.0
waveLengthX=1.0
@@ -31,7 +30,7 @@ setupTestFunction()
phaseX=0.0
phaseY=0.0
phaseZ=0.0
- sigma=0.5
+ sigma=1
# fi
}
@@ -86,11 +85,11 @@ solve()
--time-step 1 \
--time-step-order 2 \
--discrete-solver ${discreteSolver} \
- --merson-adaptivity 1.0e-7 \
+ --merson-adaptivity 1.0e-10 \
--sor-omega 1.95 \
- --gmres-restarting 20 \
- --min-iterations 20 \
- --convergence-residue 1.0e-12 \
+ --gmres-restarting 20 \
+ --min-iterations 50 \
+ --convergence-residue 1.0e-14 \
--test-function ${testFunction}\
--amplitude ${amplitude} \
--wave-length ${waveLength} \
@@ -197,4 +196,3 @@ runTest()
}
runTest
-
diff --git a/src/functors/tnlSinWaveFunction_impl.h b/src/functors/tnlSinWaveFunction_impl.h
index 60587d514a..0e64ebbdfb 100644
--- a/src/functors/tnlSinWaveFunction_impl.h
+++ b/src/functors/tnlSinWaveFunction_impl.h
@@ -137,7 +137,8 @@ getValue( const Vertex& v,
if( XDiffOrder == 0 && YDiffOrder == 2 )
return 2.0 * M_PI * y * y / ( this->waveLength * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength ) - 4.0 * M_PI * M_PI * y * y / ( this->waveLength * this->waveLength * ( x * x + y * y ) ) * this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
if( XDiffOrder == 1 && YDiffOrder == 1 )
- return -4.0 * M_PI * M_PI * x * y / ( this->waveLength * this->waveLength * x * x + y * y ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
+ return -4.0 * M_PI * M_PI * x * y / ( this->waveLength * this->waveLength * (x * x + y * y ) )* this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength )
+ - 2.0 * M_PI * this->amplitude * x * y * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength ) / ( this->waveLength * sqrt( (x * x + y * y ) * (x * x + y * y ) * (x * x + y * y ) ) );
return 0.0;
}
@@ -172,11 +173,14 @@ getValue( const Vertex& v,
if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 2 )
return 2.0 * M_PI * ( x * x + y * y ) / ( this->waveLength * sqrt( x * x + y * y + z * z ) * sqrt( x * x + y * y + z * z ) * sqrt( x * x + y * y + z * z ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength ) - 4.0 * M_PI * M_PI * z * z / ( this->waveLength * this->waveLength * ( x * x + y * y + z * z ) ) * this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength );
if( XDiffOrder == 1 && YDiffOrder == 1 && ZDiffOrder == 0 )
- return -4.0 * M_PI * M_PI * x * y / ( this->waveLength * this->waveLength * x * x + y * y + z * z ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength );
+ return -4.0 * M_PI * M_PI * x * y / ( this->waveLength * this->waveLength * (x * x + y * y + z * z ) )* this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength )
+ - 2.0 * M_PI * this->amplitude * x * y * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength ) / ( this->waveLength * sqrt( (x * x + y * y + z * z ) * (x * x + y * y + z * z ) * (x * x + y * y + z * z ) ) );
if( XDiffOrder == 1 && YDiffOrder == 0 && ZDiffOrder == 1 )
- return -4.0 * M_PI * M_PI * x * z / ( this->waveLength * this->waveLength * x * x + y * y + z * z ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength );
+ return -4.0 * M_PI * M_PI * x * z / ( this->waveLength * this->waveLength * (x * x + y * y + z * z ) )* this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength )
+ - 2.0 * M_PI * this->amplitude * x * z * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength ) / ( this->waveLength * sqrt( (x * x + y * y + z * z ) * (x * x + y * y + z * z ) * (x * x + y * y + z * z ) ) );
if( XDiffOrder == 0 && YDiffOrder == 1 && ZDiffOrder == 1 )
- return -4.0 * M_PI * M_PI * y * z / ( this->waveLength * this->waveLength * x * x + y * y + z * z ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength );
+ return -4.0 * M_PI * M_PI * z * y / ( this->waveLength * this->waveLength * (x * x + y * y + z * z ) )* this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength )
+ - 2.0 * M_PI * this->amplitude * z * y * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength ) / ( this->waveLength * sqrt( (x * x + y * y + z * z ) * (x * x + y * y + z * z ) * (x * x + y * y + z * z ) ) );
return 0.0;
}
diff --git a/src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator_impl.h b/src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator_impl.h
index 730adff00a..bc2f8da10d 100644
--- a/src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator_impl.h
+++ b/src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator_impl.h
@@ -290,12 +290,12 @@ updateLinearSystem( const RealType& time,
mesh.getHxSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, -1, 0, 0 )
+ mesh.getHySquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, 0, -1, 0 ) +
mesh.getHzSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, 0, 0, -1 ) );
- const RealType eCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time, 1, 0, 0 ) * mesh.getHxSquareInverse() /
- operatorQ.getValue(mesh, index, coordinates, u, time );
- const RealType fCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time, 0, 1, 0 ) * mesh.getHySquareInverse() /
- operatorQ.getValue(mesh, index, coordinates, u, time );
- const RealType gCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time, 0, 0, 1 ) * mesh.getHzSquareInverse() /
- operatorQ.getValue(mesh, index, coordinates, u, time );
+ const RealType eCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 1, 0, 0 );
+ const RealType fCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 0, 1, 0 );
+ const RealType gCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time * mesh.getHzSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 0, 0, 1 );
matrixRow.setElement( 0, mesh.template getCellNextToCell< 0,0,-1 >( index ), aCoef );
matrixRow.setElement( 1, mesh.template getCellNextToCell< 0,-1,0 >( index ), bCoef );
matrixRow.setElement( 2, mesh.template getCellNextToCell< -1,0,0 >( index ), cCoef );
diff --git a/src/operators/operator-Q/tnlFiniteVolumeOperatorQ_impl.h b/src/operators/operator-Q/tnlFiniteVolumeOperatorQ_impl.h
index dbc4133329..4021fd9c2a 100644
--- a/src/operators/operator-Q/tnlFiniteVolumeOperatorQ_impl.h
+++ b/src/operators/operator-Q/tnlFiniteVolumeOperatorQ_impl.h
@@ -497,19 +497,19 @@ boundaryDerivative( const MeshType& mesh,
if ( ( dx == -1 ) && ( dy == 0 ) && ( dz == 0 ) )
return mesh.getHxInverse() * ( u[ cellIndex ] - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] );
if ( ( dx == 0 ) && ( dy == 1 ) && ( dz == 0 ) )
- return mesh.getHxInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
+ return mesh.getHxInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< 1,1,0 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< -1,1,0 >( cellIndex ) ] );
if ( ( dx == 0 ) && ( dy == -1 ) && ( dz == 0 ) )
- return mesh.getHxInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
+ return mesh.getHxInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< 1,-1,0 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< -1,-1,0 >( cellIndex ) ] );
if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == 1 ) )
- return mesh.getHxInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
+ return mesh.getHxInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< 1,0,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< -1,0,1 >( cellIndex ) ] );
if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == -1 ) )
- return mesh.getHxInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
+ return mesh.getHxInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< 1,0,-1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< -1,0,-1 >( cellIndex ) ] );
}
@@ -520,19 +520,19 @@ boundaryDerivative( const MeshType& mesh,
if ( ( dx == 0 ) && ( dy == -1 ) && ( dz == 0 ) )
return mesh.getHyInverse() * ( u[ cellIndex ] - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] );
if ( ( dx == 1 ) && ( dy == 0 ) && ( dz == 0 ) )
- return mesh.getHyInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
+ return mesh.getHyInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< 1,1,0 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< 1,-1,0 >( cellIndex ) ] );
if ( ( dx == -1 ) && ( dy == 0 ) && ( dz == 0 ) )
- return mesh.getHyInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
+ return mesh.getHyInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< -1,1,0 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< -1,-1,0 >( cellIndex ) ] );
if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == 1 ) )
- return mesh.getHyInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
+ return mesh.getHyInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< 0,1,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< 0,-1,1 >( cellIndex ) ] );
if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == -1 ) )
- return mesh.getHyInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
+ return mesh.getHyInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< 0,1,-1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< 0,-1,-1 >( cellIndex ) ] );
}
@@ -543,19 +543,19 @@ boundaryDerivative( const MeshType& mesh,
if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == -1 ) )
return mesh.getHzInverse() * ( u[ cellIndex ] - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] );
if ( ( dx == 1 ) && ( dy == 0 ) && ( dz == 0 ) )
- return mesh.getHzInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
+ return mesh.getHzInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< 1,0,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< 1,0,-1 >( cellIndex ) ] );
if ( ( dx == -1 ) && ( dy == 0 ) && ( dz == 0 ) )
- return mesh.getHzInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
+ return mesh.getHzInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< -1,0,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< -1,0,-1 >( cellIndex ) ] );
if ( ( dx == 0 ) && ( dy == 1 ) && ( dz == 0 ) )
- return mesh.getHzInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
+ return mesh.getHzInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< 0,1,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< 0,1,-1 >( cellIndex ) ] );
if ( ( dx == 0 ) && ( dy == -1 ) && ( dz == 0 ) )
- return mesh.getHzInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
+ return mesh.getHzInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
u[ mesh.template getCellNextToCell< 0,-1,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] -
u[ mesh.template getCellNextToCell< 0,-1,-1 >( cellIndex ) ] );
}
--
GitLab