Commit ccd42739 authored by Jakub Klinkovský's avatar Jakub Klinkovský

Merge branch 'hamilton-jacobi-rebase' into 'develop'

Hamilton jacobi rebase

See merge request !43
parents 552e90c4 e162a57a
......@@ -36,6 +36,9 @@ class DirectEikonalSolverConfig
{
config.addDelimiter( "Direct eikonal equation solver settings:" );
config.addRequiredEntry< String >( "input-file", "Input file." );
config.addEntry< String >( "distributed-grid-io-type", "Choose Distributed Grid IO Type", "MpiIO");
config.addEntryEnum< String >( "LocalCopy" );
config.addEntryEnum< String >( "MpiIO" );
};
};
......
#!/bin/bash
device="host"
dimensions="2D 3D"
#dimensions="3D"
#dimensions="2D 3D"
dimensions="2D"
sizes1D="16 32 64 128 256 512 1024 2048 4096"
#sizes1D="256"
sizes2D="16 32 64 128 256 512 1024"
#sizes2D="8"
sizes3D="16 32 64 128 256"
sizes2D="16 32 64 128 256 512 1024 2048 4096"
#sizes2D="16"
sizes3D="8 16 32 64 128 256"
#sizes3D="128 256"
testFunctions="paraboloid"
#testFunctions="sin-wave-sdf"
#testFunctions="sin-bumps-sdf"
snapshotPeriod=0.1
finalTime=1.5
solverName="tnl-direct-eikonal-solver"
#solverName="gdb --args tnl-direct-eikonal-solver-dbg --catch-exceptions no"
#
realType="double"
#mpiRun="mpirun -np 4 -oversubscribe "
mpiRun=""
## CAREFULL: If you set LocalCopy of MPI, you have to start with mpiRun even tnl-init
## This isnt problem with MpiIO.
## CAREFULL: For smoothly calculated error, you have to choose the right output function which
## is for both MpiIO, LocalCopy different.
solverName=${mpiRun}"tnl-direct-eikonal-solver"
#solverName="gdb --args "${mpiRun}"tnl-direct-eikonal-solver-dbg --catch-exceptions no --mpi-gdb-debug false"
#scale=2.0
#finalSdf="aux-0.tnl aux-1.tnl"
finalSdf="aux-final.tnl"
setupTestFunction()
{
......@@ -22,16 +36,17 @@ setupTestFunction()
# then
origin=-1.0
proportions=2.0
# origin=-1.0
# proportions=2.0
amplitude=1.0
waveLength=0.2
waveLengthX=0.2
waveLengthY=0.2
waveLengthZ=0.2
wavesNumber=3.0
wavesNumberX=0.5
wavesNumberY=2.0
waveLength=0.4
waveLengthX=0.5
waveLengthY=0.5
waveLengthZ=0.5
wavesNumber=1.25
wavesNumberX=0.5 wavesNumberY=2.0
wavesNumberZ=3.0
phase=0.1
phase=-1.5
phaseX=0.0
phaseY=0.0
phaseZ=0.0
......@@ -44,6 +59,7 @@ setupGrid()
{
dimensions=$1
gridSize=$2
#scale=$3
tnl-grid-setup --dimensions ${dimensions} \
--origin-x ${origin} \
--origin-y ${origin} \
......@@ -53,47 +69,51 @@ setupGrid()
--proportions-z ${proportions} \
--size-x ${gridSize} \
--size-y ${gridSize} \
--size-z ${gridSize}
--real-type ${realType} \
--size-z ${gridSize}
#$((2*${gridSize}))
}
setInitialCondition()
{
testFunction=$1
tnl-init --test-function ${testFunction} \
--output-file initial-u.tnl \
--amplitude ${amplitude} \
--wave-length ${waveLength} \
--wave-length-x ${waveLengthX} \
--wave-length-y ${waveLengthY} \
--wave-length-z ${waveLengthZ} \
--waves-number ${wavesNumber} \
--waves-number-x ${wavesNumberX} \
--waves-number-y ${wavesNumberY} \
--waves-number-z ${wavesNumberZ} \
--phase ${phase} \
--phase-x ${phaseX} \
--phase-y ${phaseY} \
--phase-z ${phaseZ} \
--sigma ${sigma} \
--radius ${radius}
--real-type ${realType} \
--output-file initial-u.tnl \
--amplitude ${amplitude} \
--wave-length ${waveLength} \
--wave-length-x ${waveLengthX} \
--wave-length-y ${waveLengthY} \
--wave-length-z ${waveLengthZ} \
--waves-number ${wavesNumber} \
--waves-number-x ${wavesNumberX} \
--waves-number-y ${wavesNumberY} \
--waves-number-z ${wavesNumberZ} \
--phase ${phase} \
--phase-x ${phaseX} \
--phase-y ${phaseY} \
--phase-z ${phaseZ} \
--sigma ${sigma} \
--radius ${radius}
tnl-init --test-function ${testFunction}-sdf \
--output-file exact-u.tnl \
--amplitude ${amplitude} \
--wave-length ${waveLength} \
--wave-length-x ${waveLengthX} \
--wave-length-y ${waveLengthY} \
--wave-length-z ${waveLengthZ} \
--waves-number ${wavesNumber} \
--waves-number-x ${wavesNumberX} \
--waves-number-y ${wavesNumberY} \
--waves-number-z ${wavesNumberZ} \
--phase ${phase} \
--phase-x ${phaseX} \
--phase-y ${phaseY} \
--phase-z ${phaseZ} \
--sigma ${sigma} \
--radius ${radius}
tnl-init --test-function ${testFunction}-sdf \
--real-type ${realType} \
--output-file exact-u.tnl \
--amplitude ${amplitude} \
--wave-length ${waveLength} \
--wave-length-x ${waveLengthX} \
--wave-length-y ${waveLengthY} \
--wave-length-z ${waveLengthZ} \
--waves-number ${wavesNumber} \
--waves-number-x ${wavesNumberX} \
--waves-number-y ${wavesNumberY} \
--waves-number-z ${wavesNumberZ} \
--phase ${phase} \
--phase-x ${phaseZ} \
--phase-y ${phaseZ} \
--phase-z ${phaseZ} \
--sigma ${sigma} \
--radius ${radius} \
}
......@@ -111,17 +131,22 @@ solve()
--min-iterations 20 \
--convergence-residue 1.0e-12 \
--snapshot-period ${snapshotPeriod} \
--real-type ${realType} \
--final-time ${finalTime}
}
computeError()
{
for sweep in ${finalSdf}
do
tnl-diff --mesh mesh.tnl \
--input-files aux-final.tnl exact-u.tnl \
--input-files exact-u.tnl u-00000.tnl \
--mode sequence \
--snapshot-period ${snapshotPeriod} \
--output-file errors.txt \
--write-difference yes
#aux-final.tnl \
done
}
runTest()
......
/*
* File: tnlDirectEikonalMethodBase1D_impl.h
* Author: Fencl
*
* Created on March 15, 2019
*/
#pragma once
template< typename Real,
typename Device,
typename Index >
void
tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::
initInterface( const MeshFunctionPointer& _input,
MeshFunctionPointer& _output,
InterfaceMapPointer& _interfaceMap )
{
if( std::is_same< Device, Devices::Cuda >::value )
{
#ifdef HAVE_CUDA
const MeshType& mesh = _input->getMesh();
const int cudaBlockSize( 16 );
int numBlocksX = Devices::Cuda::getNumberOfBlocks( mesh.getDimensions().x(), cudaBlockSize );
dim3 blockSize( cudaBlockSize );
dim3 gridSize( numBlocksX );
Devices::Cuda::synchronizeDevice();
CudaInitCaller<<< gridSize, blockSize >>>( _input.template getData< Device >(),
_output.template modifyData< Device >(),
_interfaceMap.template modifyData< Device >() );
cudaDeviceSynchronize();
TNL_CHECK_CUDA_DEVICE;
#endif
}
if( std::is_same< Device, Devices::Host >::value )
{
const MeshType& mesh = _input->getMesh();
typedef typename MeshType::Cell Cell;
const MeshFunctionType& input = _input.getData();
MeshFunctionType& output = _output.modifyData();
InterfaceMapType& interfaceMap = _interfaceMap.modifyData();
Cell cell( mesh );
for( cell.getCoordinates().x() = 0;
cell.getCoordinates().x() < mesh.getDimensions().x();
cell.getCoordinates().x() ++ )
{
cell.refresh();
output[ cell.getIndex() ] =
input( cell ) >= 0 ? std::numeric_limits< RealType >::max() :
-std::numeric_limits< RealType >::max();
interfaceMap[ cell.getIndex() ] = false;
}
const RealType& h = mesh.getSpaceSteps().x();
for( cell.getCoordinates().x() = 0;
cell.getCoordinates().x() < mesh.getDimensions().x() - 1;
cell.getCoordinates().x() ++ )
{
cell.refresh();
const RealType& c = input( cell );
if( ! cell.isBoundaryEntity() )
{
const auto& neighbors = cell.getNeighborEntities();
Real pom = 0;
//const IndexType& c = cell.getIndex();
const IndexType e = neighbors.template getEntityIndex< 1 >();
if( c * input[ e ] <= 0 )
{
pom = TNL::sign( c )*( h * c )/( c - input[ e ]);
if( TNL::abs( output[ cell.getIndex() ] ) > TNL::abs( pom ) )
output[ cell.getIndex() ] = pom;
pom = pom - TNL::sign( c )*h; //output[ e ] = (hx * c)/( c - input[ e ]) - hx;
if( TNL::abs( output[ e ] ) > TNL::abs( pom ) )
output[ e ] = pom;
interfaceMap[ cell.getIndex() ] = true;
interfaceMap[ e ] = true;
}
}
}
}
}
template< typename Real,
typename Device,
typename Index >
template< typename MeshEntity >
void
tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::
updateCell( MeshFunctionType& u,
const MeshEntity& cell,
const RealType v )
{
const auto& neighborEntities = cell.template getNeighborEntities< 1 >();
const MeshType& mesh = cell.getMesh();
const RealType& h = mesh.getSpaceSteps().x();
const RealType value = u( cell );
RealType a, tmp = std::numeric_limits< RealType >::max();
if( cell.getCoordinates().x() == 0 )
a = u[ neighborEntities.template getEntityIndex< 1 >() ];
else if( cell.getCoordinates().x() == mesh.getDimensions().x() - 1 )
a = u[ neighborEntities.template getEntityIndex< -1 >() ];
else
{
a = TNL::argAbsMin( u[ neighborEntities.template getEntityIndex< -1 >() ],
u[ neighborEntities.template getEntityIndex< 1 >() ] );
}
if( fabs( a ) == std::numeric_limits< RealType >::max() )
return;
tmp = a + TNL::sign( value ) * h/v;
u[ cell.getIndex() ] = argAbsMin( value, tmp );
}
template< typename Real,
typename Device,
typename Index >
__cuda_callable__
bool
tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::
updateCell( volatile Real sArray[18], int thri, const Real h, const Real v )
{
const RealType value = sArray[ thri ];
RealType a, tmp = std::numeric_limits< RealType >::max();
a = TNL::argAbsMin( sArray[ thri+1 ],
sArray[ thri-1 ] );
if( fabs( a ) == std::numeric_limits< RealType >::max() )
return false;
tmp = a + TNL::sign( value ) * h/v;
sArray[ thri ] = argAbsMin( value, tmp );
tmp = value - sArray[ thri ];
if ( fabs( tmp ) > 0.001*h )
return true;
else
return false;
}
#ifdef HAVE_CUDA
template < typename Real, typename Device, typename Index >
__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 1, Real, Device, Index > >& input,
Functions::MeshFunction< Meshes::Grid< 1, Real, Device, Index > >& output,
Functions::MeshFunction< Meshes::Grid< 1, Real, Device, Index >, 1, bool >& interfaceMap )
{
int i = threadIdx.x + blockDim.x*blockIdx.x;
const Meshes::Grid< 1, Real, Device, Index >& mesh = input.template getMesh< Devices::Cuda >();
if( i < mesh.getDimensions().x() )
{
typedef typename Meshes::Grid< 1, Real, Device, Index >::Cell Cell;
Cell cell( mesh );
cell.getCoordinates().x() = i;
cell.refresh();
const Index cind = cell.getIndex();
output[ cind ] =
input( cell ) >= 0 ? std::numeric_limits< Real >::max() :
- std::numeric_limits< Real >::max();
interfaceMap[ cind ] = false;
const Real& h = mesh.getSpaceSteps().x();
cell.refresh();
const Real& c = input( cell );
if( ! cell.isBoundaryEntity() )
{
auto neighbors = cell.getNeighborEntities();
Real pom = 0;
const Index e = neighbors.template getEntityIndex< 1 >();
const Index w = neighbors.template getEntityIndex< -1 >();
if( c * input[ e ] <= 0 )
{
pom = TNL::sign( c )*( h * c )/( c - input[ e ]);
if( TNL::abs( output[ cind ] ) > TNL::abs( pom ) )
output[ cind ] = pom;
interfaceMap[ cind ] = true;
}
if( c * input[ w ] <= 0 )
{
pom = TNL::sign( c )*( h * c )/( c - input[ w ]);
if( TNL::abs( output[ cind ] ) > TNL::abs( pom ) )
output[ cind ] = pom;
interfaceMap[ cind ] = true;
}
}
}
}
#endif
/*
* File: tnlDirectEikonalMethodBase2D_impl.h
* Author: Fencl
*
* Created on March 15, 2019
*/
#pragma once
template< typename Real,
typename Device,
typename Index >
void
tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >::
initInterface( const MeshFunctionPointer& _input,
MeshFunctionPointer& _output,
InterfaceMapPointer& _interfaceMap,
const StaticVector vecLowerOverlaps,
const StaticVector vecUpperOverlaps )
{
if( std::is_same< Device, Devices::Cuda >::value )
{
#ifdef HAVE_CUDA
const MeshType& mesh = _input->getMesh();
const int cudaBlockSize( 16 );
int numBlocksX = Devices::Cuda::getNumberOfBlocks( mesh.getDimensions().x(), cudaBlockSize );
int numBlocksY = Devices::Cuda::getNumberOfBlocks( mesh.getDimensions().y(), cudaBlockSize );
dim3 blockSize( cudaBlockSize, cudaBlockSize );
dim3 gridSize( numBlocksX, numBlocksY );
Devices::Cuda::synchronizeDevice();
CudaInitCaller<<< gridSize, blockSize >>>( _input.template getData< Device >(),
_output.template modifyData< Device >(),
_interfaceMap.template modifyData< Device >(),
vecLowerOverlaps, vecUpperOverlaps);
cudaDeviceSynchronize();
TNL_CHECK_CUDA_DEVICE;
#endif
}
if( std::is_same< Device, Devices::Host >::value )
{
MeshFunctionType input = _input.getData();
MeshFunctionType& output = _output.modifyData();
InterfaceMapType& interfaceMap = _interfaceMap.modifyData();
const MeshType& mesh = input.getMesh();
typedef typename MeshType::Cell Cell;
Cell cell( mesh );
for( cell.getCoordinates().y() = 0;
cell.getCoordinates().y() < mesh.getDimensions().y();
cell.getCoordinates().y() ++ )
for( cell.getCoordinates().x() = 0;
cell.getCoordinates().x() < mesh.getDimensions().x();
cell.getCoordinates().x() ++ )
{
cell.refresh();
output[ cell.getIndex() ] =
input( cell ) >= 0 ? std::numeric_limits< RealType >::max() :
- std::numeric_limits< RealType >::max();
interfaceMap[ cell.getIndex() ] = false;
}
const RealType& hx = mesh.getSpaceSteps().x();
const RealType& hy = mesh.getSpaceSteps().y();
for( cell.getCoordinates().y() = 0 + vecLowerOverlaps[1];
cell.getCoordinates().y() < mesh.getDimensions().y() - vecUpperOverlaps[1];
cell.getCoordinates().y() ++ )
for( cell.getCoordinates().x() = 0 + vecLowerOverlaps[0];
cell.getCoordinates().x() < mesh.getDimensions().x() - vecUpperOverlaps[0];
cell.getCoordinates().x() ++ )
{
cell.refresh();
const RealType& c = input( cell );
if( ! cell.isBoundaryEntity() )
{
auto neighbors = cell.getNeighborEntities();
Real pom = 0;
const IndexType e = neighbors.template getEntityIndex< 1, 0 >();
const IndexType n = neighbors.template getEntityIndex< 0, 1 >();
if( c * input[ n ] <= 0 )
{
pom = TNL::sign( c )*( hy * c )/( c - input[ n ]);
if( TNL::abs( output[ cell.getIndex() ] ) > TNL::abs( pom ) )
output[ cell.getIndex() ] = pom;
pom = pom - TNL::sign( c )*hy;
if( TNL::abs( output[ n ] ) > TNL::abs( pom ) )
output[ n ] = pom; //( hy * c )/( c - input[ n ]) - hy;
interfaceMap[ cell.getIndex() ] = true;
interfaceMap[ n ] = true;
}
if( c * input[ e ] <= 0 )
{
pom = TNL::sign( c )*( hx * c )/( c - input[ e ]);
if( TNL::abs( output[ cell.getIndex() ] ) > TNL::abs( pom ) )
output[ cell.getIndex() ] = pom;
pom = pom - TNL::sign( c )*hx; //output[ e ] = (hx * c)/( c - input[ e ]) - hx;
if( TNL::abs( output[ e ] ) > TNL::abs( pom ) )
output[ e ] = pom;
interfaceMap[ cell.getIndex() ] = true;
interfaceMap[ e ] = true;
}
}
}
}
}
template< typename Real,
typename Device,
typename Index >
template< typename MeshEntity >
__cuda_callable__
bool
tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >::
updateCell( MeshFunctionType& u,
const MeshEntity& cell,
const RealType v)
{
const auto& neighborEntities = cell.template getNeighborEntities< 2 >();
const MeshType& mesh = cell.getMesh();
const RealType& hx = mesh.getSpaceSteps().x();
const RealType& hy = mesh.getSpaceSteps().y();
const RealType value = u( cell );
RealType a, b, tmp = std::numeric_limits< RealType >::max();
if( cell.getCoordinates().x() == 0 )
a = u[ neighborEntities.template getEntityIndex< 1, 0 >() ];
else if( cell.getCoordinates().x() == mesh.getDimensions().x() - 1 )
a = u[ neighborEntities.template getEntityIndex< -1, 0 >() ];
else
{
a = TNL::argAbsMin( u[ neighborEntities.template getEntityIndex< -1, 0 >() ],
u[ neighborEntities.template getEntityIndex< 1, 0 >() ] );
}
if( cell.getCoordinates().y() == 0 )
b = u[ neighborEntities.template getEntityIndex< 0, 1 >()];
else if( cell.getCoordinates().y() == mesh.getDimensions().y() - 1 )
b = u[ neighborEntities.template getEntityIndex< 0, -1 >() ];
else
{
b = TNL::argAbsMin( u[ neighborEntities.template getEntityIndex< 0, -1 >() ],
u[ neighborEntities.template getEntityIndex< 0, 1 >() ] );
}
if( fabs( a ) == std::numeric_limits< RealType >::max() &&
fabs( b ) == std::numeric_limits< RealType >::max() )
return false;
RealType pom[6] = { a, b, std::numeric_limits< RealType >::max(), hx, hy, 0.0 };
tmp = getNewValue( pom , value, v );
u[ cell.getIndex() ] = tmp;
tmp = value - u[ cell.getIndex() ];
if ( fabs( tmp ) > 0.001*hx )
return true;
else
return false;
}
template< typename Real,
typename Device,
typename Index >
template< int sizeSArray >
__cuda_callable__
bool
tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >::
updateCell( volatile Real *sArray, int thri, int thrj, const Real hx, const Real hy,
const Real v )
{
const RealType value = sArray[ thrj * sizeSArray + thri ];
RealType a, b, tmp = std::numeric_limits< RealType >::max();
b = TNL::argAbsMin( sArray[ (thrj+1) * sizeSArray + thri ],
sArray[ (thrj-1) * sizeSArray + thri ] );
a = TNL::argAbsMin( sArray[ thrj * sizeSArray + thri+1 ],
sArray[ thrj * sizeSArray + thri-1 ] );
if( fabs( a ) == std::numeric_limits< RealType >::max() &&
fabs( b ) == std::numeric_limits< RealType >::max() )
return false;
RealType pom[6] = { a, b, std::numeric_limits< RealType >::max(), (RealType)hx, (RealType)hy, 0.0 };
tmp = getNewValue( pom , value, v );
sArray[ thrj * sizeSArray + thri ] = tmp;
tmp = value - sArray[ thrj * sizeSArray + thri ];
if ( fabs( tmp ) > 0.001*hx )
return true;
else
return false;
}
template< typename Real,
typename Device,
typename Index >
__cuda_callable__
Real
tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >::
getNewValue( RealType valuesAndSteps[], const RealType originalValue, const RealType v )
{
RealType newValue = std::numeric_limits< RealType >::max();
sortMinims( valuesAndSteps );
// calculation of real value taken from ZHAO
newValue = valuesAndSteps[ 0 ] + TNL::sign( originalValue ) * valuesAndSteps[ 3 ]/v;
if( fabs( newValue ) < fabs( valuesAndSteps[ 1 ] ) )
{
newValue = argAbsMin( originalValue, newValue );
}
else
{
newValue = ( valuesAndSteps[ 3 ] * valuesAndSteps[ 3 ] * valuesAndSteps[ 1 ] +
valuesAndSteps[ 4 ] * valuesAndSteps[ 4 ] * valuesAndSteps[ 0 ] +
TNL::sign( originalValue ) * valuesAndSteps[ 3 ] * valuesAndSteps[ 4 ] *
TNL::sqrt( ( valuesAndSteps[ 3 ] * valuesAndSteps[ 3 ] + valuesAndSteps[ 4 ] * valuesAndSteps[ 4 ] )/( v * v ) -
( valuesAndSteps[ 1 ] - valuesAndSteps[ 0 ] ) *
( valuesAndSteps[ 1 ] - valuesAndSteps[ 0 ] ) ) )/
( valuesAndSteps[ 3 ] * valuesAndSteps[ 3 ] + valuesAndSteps[ 4 ] * valuesAndSteps[ 4 ] );
newValue = argAbsMin( originalValue, newValue );
}
return newValue;
}
template < typename T1 >
__cuda_callable__ void sortMinims( T1 pom[] )
{
T1 tmp[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
if( fabs(pom[0]) <= fabs(pom[1]) && fabs(pom[1]) <= fabs(pom[2])){
tmp[0] = pom[0]; tmp[1] = pom[1]; tmp[2] = pom[2];
tmp[3] = pom[3]; tmp[4] = pom[4]; tmp[5] = pom[5];
}
else if( fabs(pom[0]) <= fabs(pom[2]) && fabs(pom[2]) <= fabs(pom[1]) ){
tmp[0] = pom[0]; tmp[1] = pom[2]; tmp[2] = pom[1];
tmp[3] = pom[3]; tmp<