Commit 48f3bfa7 authored by Quentin Khan's avatar Quentin Khan
Browse files

Merge branch 'master' into adaptive

parents a70cdd5d 4941a7b8
......@@ -115,12 +115,12 @@ public:
}else{
if(type==SOURCE){
for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint<FReal>(&XYZ[3*idPart]),FParticleTypeSource,idPart);
octree->insert(FPoint<FReal>(&XYZ[3*idPart]),FParticleType::FParticleTypeSource,idPart);
}
FScalFMMEngine<FReal>::nbPart += NbPositions;
}else{
for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint<FReal>(&XYZ[3*idPart]),FParticleTypeTarget,idPart);
octree->insert(FPoint<FReal>(&XYZ[3*idPart]),FParticleType::FParticleTypeTarget,idPart);
}
FScalFMMEngine<FReal>::nbPart += NbPositions;
}
......@@ -138,12 +138,12 @@ public:
}else{
if(type==SOURCE){
for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),FParticleTypeSource,idPart);
octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),FParticleType::FParticleTypeSource,idPart);
}
FScalFMMEngine<FReal>::nbPart += NbPositions;
}else{
for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),FParticleTypeTarget,idPart);
octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),FParticleType::FParticleTypeTarget,idPart);
}
FScalFMMEngine<FReal>::nbPart += NbPositions;
}
......
......@@ -380,12 +380,12 @@ public:
}else{
if(type==SOURCE){
for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),FParticleTypeSource,idPart);
octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),FParticleType::FParticleTypeSource,idPart);
}
FScalFMMEngine<FReal>::nbPart += NbPositions;
}else{
for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),FParticleTypeTarget,idPart);
octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),FParticleType::FParticleTypeTarget,idPart);
}
FScalFMMEngine<FReal>::nbPart += NbPositions;
}
......@@ -403,12 +403,12 @@ public:
}else{
if(type==SOURCE){
for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint<FReal>(&XYZ[3*idPart]),FParticleTypeSource,idPart);
octree->insert(FPoint<FReal>(&XYZ[3*idPart]),FParticleType::FParticleTypeSource,idPart);
}
FScalFMMEngine<FReal>::nbPart += NbPositions;
}else{
for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint<FReal>(&XYZ[3*idPart]),FParticleTypeTarget,idPart);
octree->insert(FPoint<FReal>(&XYZ[3*idPart]),FParticleType::FParticleTypeTarget,idPart);
}
FScalFMMEngine<FReal>::nbPart += NbPositions;
}
......
......@@ -321,6 +321,13 @@ if (MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/CMakeModules/morse/
##############################################################################
#
if( SCALFMM_USE_BLAS )
# include(FortranCInterface)
# # Define a Fortran interface file (FCMangle.hpp)
# FortranCInterface_HEADER( ${CMAKE_CURRENT_SOURCE_DIR}/Src/FCMangle.hpp
# MACRO_NAMESPACE "PM_"
# SYMBOL_NAMESPACE "PM_"
# SYMBOLS init testPPM:init)
set(SCALFMM_BLAS_ADD_ ON) # A enlever lorsque MORSE detectera le mangling fortan
message(STATUS "CMAKE_CXX_COMPILER_ID STREQUAL ${CMAKE_CXX_COMPILER_ID}")
option( SCALFMM_USE_MKL_AS_BLAS "Set to ON to use MKL CBLAS" OFF )
......
8 4
20000
0.5 0.5 0.5 0.5
20000 1 0.5 0.5 0.5
0.840188 0.394383 0.783099 0.01 1
0.911647 0.197551 0.335223 0.01 1
0.277775 0.55397 0.477397 0.01 1
......@@ -20000,4 +19998,4 @@
0.00448784 0.00539908 0.182474 0.01 0
0.0237434 0.139661 0.412617 0.01 1
0.514349 0.627817 0.0209046 0.01 1
0.56572 0.990817 0.904442 0.01 0
\ No newline at end of file
0.56572 0.990817 0.904442 0.01 0
......@@ -107,9 +107,9 @@ If you consider the
\subsection{Plummer Model}
This is a hard test case in astrophysics problem, and it models a globular cluster of stars, which is highly non uniform. It is called the plummer distribution. To construct such distribution, first we construct a uniform points distribution on the unit sphere. Second, the radius is chosen according to the plummer distribution (double power law in astrophysics). We consider $u$ a random number between 0 and 1, then the associated radius is given by
\begin{equation*}
r = \sqrt{\frac{u^{2/3}}{u^{2/3}-1}}
r = 1.0/\sqrt{u^{-2/3}-1},
\end{equation*}
and the total mass is one. Then, $m_i = \frac{1}{Npt}$.
\begin{figure}[h]
\centering
\begin{minipage}{0.45\textwidth}%
......@@ -140,6 +140,8 @@ The corresponding potential is
\begin{equation}
\Phi_P(r) = - \frac{G M}{\sqrt{r^2+a^2}}
\end{equation}
In N-body units, $G = M = 1$ and $a = 3\pi/16 \sim 0.589$
\subsection{Diagonal Model}
%, shape end size=.5cm},decoration={shape start size=.5cm, shape end size=.125cm
......
......@@ -97,7 +97,7 @@ namespace Param {
const FParameterNames Size
= {{"-size"}, "Size of the geometry a:b:c - default values are 1.0:1.0:2.0"};
const FParameterNames Radius
= {{"-radius"}, "used to specified the radius of the sphere an dthe plummer distribution or R - default value for R is 2.0"};
= {{"-radius"}, "used to specified the radius of the sphere and the plummer distribution or R - default value for R is 2.0"};
const FParameterNames Charge
= {{"-charge"}, "generate physical values between -1 and 1 otherwise generate between 0 and 1"};
const FParameterNames ZM
......@@ -250,6 +250,10 @@ int main(int argc, char ** argv){
std::cout << "End of writing" <<std::endl;
// Generate file for visualization
// if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputVisuFile.options)){
// std::string outfilename(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.vtp"));
// driverExportData(outfilename, particles , NbPoints,loader.getNbRecordPerline() );
// }
if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputVisuFile.options)) {
std::string visufile(FParameters::getStr(argc, argv, FParameterDefinitions::OutputVisuFile.options, "output.vtp"));
driverExportData(visufile, particles , NbPoints);
......
......@@ -85,7 +85,7 @@ public:
const MortonIndex particuleIndex = tree->getMortonFromPosition(currentPart);
if(particuleIndex != currentMortonIndex){
//Need to move this one
interface->removeFromLeafAndKeep(particles,currentPart,idxPart,FParticleTypeSource);
interface->removeFromLeafAndKeep(particles,currentPart,idxPart,FParticleType::FParticleTypeSource);
}
else{
//Need to increment idx;
......@@ -102,7 +102,7 @@ public:
const MortonIndex particuleIndex = tree->getMortonFromPosition(currentPart);
if(particuleIndex != currentMortonIndex){
//Need to move this one
interface->removeFromLeafAndKeep(particleTargets,currentPart,idxPart, FParticleTypeTarget);
interface->removeFromLeafAndKeep(particleTargets,currentPart,idxPart, FParticleType::FParticleTypeTarget);
}
else{
//Need to increment idx;
......
......@@ -33,11 +33,11 @@ public:
for(int idxAttr = 0 ; idxAttr < ContainerClass::NbAttributes ; ++idxAttr){
particleValues[idxAttr] = lf->getAttribute(idxAttr)[idxPart];
}
if(type == FParticleTypeTarget){
toStoreRemovedTargetParts.push(particlePos,FParticleTypeTarget,lf->getIndexes()[idxPart],particleValues);
if(type == FParticleType::FParticleTypeTarget){
toStoreRemovedTargetParts.push(particlePos,FParticleType::FParticleTypeTarget,lf->getIndexes()[idxPart],particleValues);
}
else{
toStoreRemovedSourceParts.push(particlePos,FParticleTypeSource,lf->getIndexes()[idxPart],particleValues);
toStoreRemovedSourceParts.push(particlePos,FParticleType::FParticleTypeSource,lf->getIndexes()[idxPart],particleValues);
}
lf->removeParticles(&idxPart,1);
}
......@@ -53,7 +53,7 @@ public:
const FPoint<FReal> particlePos(toStoreRemovedSourceParts.getPositions()[0][idxToInsert],
toStoreRemovedSourceParts.getPositions()[1][idxToInsert],
toStoreRemovedSourceParts.getPositions()[2][idxToInsert]);
tree->insert(particlePos, FParticleTypeSource, toStoreRemovedSourceParts.getIndexes()[idxToInsert], particleValues);
tree->insert(particlePos, FParticleType::FParticleTypeSource, toStoreRemovedSourceParts.getIndexes()[idxToInsert], particleValues);
}
for(FSize idxToInsert = 0; idxToInsert<toStoreRemovedTargetParts.getNbParticles() ; ++idxToInsert){
......@@ -64,7 +64,7 @@ public:
toStoreRemovedTargetParts.getPositions()[1][idxToInsert],
toStoreRemovedTargetParts.getPositions()[2][idxToInsert]);
tree->insert(particlePos, FParticleTypeTarget, toStoreRemovedTargetParts.getIndexes()[idxToInsert], particleValues);
tree->insert(particlePos, FParticleType::FParticleTypeTarget, toStoreRemovedTargetParts.getIndexes()[idxToInsert], particleValues);
}
toStoreRemovedSourceParts.clear();
......
......@@ -19,7 +19,7 @@
/**
* @brief The FParticleType enum is to make a difference between Target and Source (Tsm)
*/
enum FParticleType {
enum class FParticleType {
FParticleTypeSource = 0,
FParticleTypeTarget = 1
};
......
......@@ -51,8 +51,8 @@ public:
*/
template<typename... Args>
void push(const FPoint<FReal>& inParticlePosition, const FParticleType type, Args ... args){
if(type == FParticleTypeTarget) targets.push(inParticlePosition, FParticleTypeTarget, args...);
else sources.push(inParticlePosition, FParticleTypeSource, args...);
if(type == FParticleType::FParticleTypeTarget) targets.push(inParticlePosition, FParticleType::FParticleTypeTarget, args...);
else sources.push(inParticlePosition, FParticleType::FParticleTypeSource, args...);
}
/**
......
......@@ -131,8 +131,8 @@ public:
inParticlePositions->setPosition(x,y,z);
*inPhysicalValue = data;
if(isTarget) (*particleType) = FParticleTypeTarget;
else (*particleType) = FParticleTypeSource;
if(isTarget) (*particleType) = FParticleType::FParticleTypeTarget;
else (*particleType) = FParticleType::FParticleTypeSource;
}
};
......
......@@ -243,12 +243,12 @@ FReal plummerDist(FSize& cpt, const FReal &R) {
}
}
/**
* \brief Build N points following the Plummer distribution
*
* First we construct N points uniformly distributed on the unit sphere. Then
* the radius in construct according to the Plummr distribution.
* the radius in construct according to the Plummer distribution for
* a constant mass of 1/N
*
* \tparam FReal Floating point type
*
......@@ -259,6 +259,32 @@ FReal plummerDist(FSize& cpt, const FReal &R) {
template <class FReal>
void unifRandomPlummer(const FSize N, const FReal R, FReal * points) {
unifRandomPointsOnSphere<FReal>(N, 1, points);
FReal mc = 1.0/static_cast<FReal>(N);
for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4) {
FReal m = getRandom<FReal>();
FReal r = FMath::Sqrt( 1.0/(FMath::pow(m, -2.0/3.0) - 1.0)) ;
points[j] *= r;
points[j+1] *= r;
points[j+2] *= r;
points[j+3] = mc; // the mass
}
}
/**
* \brief Build N points following the Plummer like distribution
*
* First we construct N points uniformly distributed on the unit sphere. Then
* the radius in construct according to the Plummer like distribution.
*
* \tparam FReal Floating point type
*
* \param N the number of points following the Plummer distribution
* \param R the radius of the sphere that contains all the points
* \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
*/
template <class FReal>
void unifRandomPlummerLike(const FSize N, const FReal R, FReal * points) {
FReal a = 1.0 ;
unifRandomPointsOnSphere<FReal>(N, 1, points);
FReal r, rm = 0.0;
FSize cpt = 0;
for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4) {
......
......@@ -124,8 +124,8 @@ public:
void fillParticle(FPoint<FReal>*const inParticlePositions, FParticleType*const isTarget){
FRandomLoader<FReal>::fillParticle(inParticlePositions);
if(FRandomLoader<FReal>::getRandom() > 0.5 ) (*isTarget) = FParticleTypeTarget;
else (*isTarget) = FParticleTypeSource;
if(FRandomLoader<FReal>::getRandom() > 0.5 ) (*isTarget) = FParticleType::FParticleTypeTarget;
else (*isTarget) = FParticleType::FParticleTypeSource;
}
};
......
......@@ -28,7 +28,10 @@
#cmakedefine SCALFMM_USE_BLAS
#cmakedefine SCALFMM_USE_MKL_AS_BLAS
// Fortran Mangling
#cmakedefine SCALFMM_BLAS_ADD_
#cmakedefine SCALFMM_BLAS_UPCASE
#cmakedefine SCALFMM_BLAS_NOCHANGE
////////////////////////////////////////////////////////
// FFT
///////////////////////////////////////////////////////
......
This diff is collapsed.
/*
* FFortranMangling.hpp
*
* Created on: 6 juin 2016
* Author: coulaud
*/
#ifndef SRC_UTILS_FFORTRANMANGLING_HPP_
#define SRC_UTILS_FFORTRANMANGLING_HPP_
#include "ScalFmmConfig.h"
#ifdef SCALFMM_BLAS_ADD_
/* Mangling for Fortran subroutine symbols with underscores. */
#define FortranName(name,NAME) name##_
#elif defined(SCALFMM_BLAS_UPCASE)
/* Mangling for Fortran subroutine symbols in uppercase and without underscores */
#define FortranName(name,NAME) NAME
#elif defined(SCALFMM_BLAS_NOCHANGE)
/* Mangling for Fortran subroutine symbols without no change. */
#define FortranName(name,NAME) name
#else
#error("Fortran MANGLING NOT DEFINED")
#endif
// blas 1
#define Fddot FortranName(ddot,DDOT)
#define Fdscal FortranName(dscal,DSCAL)
#define Fdcopy FortranName(dcopy,DCOPY)
#define Fdaxpy FortranName(daxpy,DAXPY)
#define Fsdot FortranName(sdot,SDOT)
#define Fsscal FortranName(sscal,SSCAL)
#define Fscopy FortranName(scopy,SCOPY)
#define Fsaxpy FortranName(saxpy,SAXPY)
#define Fcscal FortranName(cscal,CSCAL)
#define Fccopy FortranName(ccopy,CCOPY)
#define Fcaxpy FortranName(caxpy,CAXPY)
#define Fzscal FortranName(zscal,ZSCAL)
#define Fzcopy FortranName(zcopy,ZCOPY)
#define Fzaxpy FortranName(zaxpy,ZAXPY)
// blas 2
#define Fdgemv FortranName(dgemv,DGEMV)
#define Fsgemv FortranName(sgemv,SGEMV)
#define Fcgemv FortranName(cgemv,CGEMV)
#define Fzgemv FortranName(zgemv,ZGEMV)
// blas 3
#define Fdgemm FortranName(dgemm,DGEMM)
#define Fsgemm FortranName(sgemm,SGEMM)
#define Fcgemm FortranName(cgemm,CGEMM)
#define Fzgemm FortranName(zgemm,ZGEMM)
// lapack
#define Fdgesvd FortranName(dgesvd,DGESVD)
#define Fdgeqrf FortranName(dgeqrf,DGEQRF)
#define Fdgeqp3 FortranName(dgeqp3,DGEQP3)
#define Fdorgqr FortranName(dorgqr,DORGQR)
#define Fdormqr FortranName(dormqr,DORMQR)
#define Fdpotrf FortranName(dpotrf,DPOTRF)
#define Fsgesvd FortranName(sgesvd,SGESVD)
#define Fsgeqrf FortranName(sgeqrf,SGEQRF)
#define Fsgeqp3 FortranName(sgeqp3,SGEQP3)
#define Fsorgqr FortranName(sorgqr,SORGQR)
#define Fsormqr FortranName(sormqr,SORMQR)
#define Fspotrf FortranName(spotrf,SPOTRF)
#define Fcgesvd FortranName(cgesvd,CGESVD)
#define Fcgeqrf FortranName(cgeqrf,CGEQRF)
#define Fcgeqp3 FortranName(cgeqp3,CGEQP3)
#define Fcorgqr FortranName(corgqr,CORGQR)
#define Fcormqr FortranName(cormqr,CORMQR)
#define Fcpotrf FortranName(cpotrf,CPOTRF)
#define Fzgesvd FortranName(zgesvd,ZGESVD)
#define Fzgeqrf FortranName(zgeqrf,ZGEQRF)
#define Fzgeqp3 FortranName(zgeqp3,ZGEQP3)
#define Fzorgqr FortranName(zorgqr,ZORGQR)
#define Fzormqr FortranName(zormqr,ZORMQR)
#define Fzpotrf FortranName(zpotrf,ZPOTRF)
#endif /* SRC_UTILS_FFORTRANMANGLING_HPP_ */
......@@ -98,7 +98,7 @@ int main(int argc, char ** argv){
(BoxWidth*FReal(drand48())) + (BoxCenter-(BoxWidth/2)),
(BoxWidth*FReal(drand48())) + (BoxCenter-(BoxWidth/2)),
(BoxWidth*FReal(drand48())) + (BoxCenter-(BoxWidth/2)));
tree.insert(particleToFill,FParticleTypeSource,idxPart);
tree.insert(particleToFill,FParticleType::FParticleTypeSource,idxPart);
}
for(FSize idxPart = 0 ; idxPart < NbPart_Target; ++idxPart){
......@@ -106,7 +106,7 @@ int main(int argc, char ** argv){
(BoxWidth*FReal(drand48())) + (BoxCenter-(BoxWidth/2)),
(BoxWidth*FReal(drand48())) + (BoxCenter-(BoxWidth/2)),
(BoxWidth*FReal(drand48())) + (BoxCenter-(BoxWidth/2)));
tree.insert(particleToFill,FParticleTypeTarget,idxPart);
tree.insert(particleToFill,FParticleType::FParticleTypeTarget,idxPart);
}
}
......
......@@ -69,7 +69,7 @@ class TestChebyshevDirectTsm : public FUTester<TestChebyshevDirectTsm> {
FPoint<FReal> position;
loader.fillParticle(&position);
// put in tree
tree.insert(position, FParticleTypeTarget, idxPart, physicalValue);
tree.insert(position, FParticleType::FParticleTypeTarget, idxPart, physicalValue);
// get copy
particlesTargets[idxPart].setPosition(position);
*(particlesTargets[idxPart].setPhysicalValue()) = physicalValue;
......@@ -84,7 +84,7 @@ class TestChebyshevDirectTsm : public FUTester<TestChebyshevDirectTsm> {
FPoint<FReal> position;
loader.fillParticle(&position);
// put in tree
tree.insert(position, FParticleTypeSource, idxPart, physicalValue);
tree.insert(position, FParticleType::FParticleTypeSource, idxPart, physicalValue);
// get copy
particlesSources[idxPart].setPosition(position);
*(particlesSources[idxPart].setPhysicalValue()) = physicalValue;
......
......@@ -68,7 +68,7 @@ class TestRotationDirectTsm : public FUTester<TestRotationDirectTsm> {
FPoint<FReal> position;
loader.fillParticle(&position);
// put in tree
tree.insert(position, FParticleTypeTarget, idxPart, physicalValue);
tree.insert(position, FParticleType::FParticleTypeTarget, idxPart, physicalValue);
// get copy
particlesTargets[idxPart].setPosition(position);
*(particlesTargets[idxPart].setPhysicalValue()) = physicalValue;
......@@ -83,7 +83,7 @@ class TestRotationDirectTsm : public FUTester<TestRotationDirectTsm> {
FPoint<FReal> position;
loader.fillParticle(&position);
// put in tree
tree.insert(position, FParticleTypeSource, idxPart, physicalValue);
tree.insert(position, FParticleType::FParticleTypeSource, idxPart, physicalValue);
// get copy
particlesSources[idxPart].setPosition(position);
*(particlesSources[idxPart].setPhysicalValue()) = physicalValue;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment