Commit 6f005bcf authored by Ondrej Szekely's avatar Ondrej Szekely
Browse files

THE LAST FIXING!!!!

parent e3943328

examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test

+12 −14
Original line number Diff line number Diff line 
#!/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 \

                 --min-iterations 50 \

                 --convergence-residue 1.0e-14 \

                 --test-function ${testFunction}\

                 --amplitude ${amplitude} \

                 --wave-length ${waveLength} \
@@ -197,4 +196,3 @@ runTest() 
}



runTest

src/functors/tnlSinWaveFunction_impl.h

+8 −4
Original line number Diff line number Diff line
@@ -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;

}

src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator_impl.h

+6 −6
Original line number Diff line number Diff line
@@ -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 );

src/operators/operator-Q/tnlFiniteVolumeOperatorQ_impl.h

+12 −12
Original line number Diff line number Diff line
@@ -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 ) ] );

    }