diff --git a/CMakeLists.txt b/CMakeLists.txt
index 77936e8f3bc4c8c7beaaab9567e8ab1daae9ba78..6da6b982700a3c38b5a8f96dec09892f7fe4e2de 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -194,6 +194,15 @@ endif()
if(WITH_EXTRAS)
# Arnoldi solver needs zgeev, which means lapack
compile_with(LAPACK)
+ foreach(_file ${LAPACK_LIBRARY})
+ get_filename_component(_name ${_file} DIRECTORY)
+ list(FIND dirlist ${_name} isin)
+ if(isin EQUAL -1)
+ list(APPEND dirlist ${_name})
+ endif()
+ endforeach()
+ set(EXTRASLIB ${dirlist} CACHE PATH "extras libraries dir")
+
endif()
# ========= Check which opencl devices are available on the system =========
@@ -368,10 +377,10 @@ if(USE_FORTRAN)
endif()
if(USE_CXX OR USE_FORTRAN)
- #add_subdirectory(src)
- #get_directory_property(FORTRAN_INCLUDE_DIRS
- #DIRECTORY ${CMAKE_SOURCE_DIR}/src
- #DEFINITION FORTRAN_INCLUDE_DIRS)
+ add_subdirectory(src)
+ get_directory_property(FORTRAN_INCLUDE_DIRS
+ DIRECTORY ${CMAKE_SOURCE_DIR}/src
+ DEFINITION FORTRAN_INCLUDE_DIRS)
endif()
if(USE_CXX)
diff --git a/cmake/fortran_utils.cmake b/cmake/fortran_utils.cmake
index b98468b629740e6c50cf5100c29f04efdd03e259..6805925b6628dc73952d090b5ae228e07d51ee0b 100644
--- a/cmake/fortran_utils.cmake
+++ b/cmake/fortran_utils.cmake
@@ -31,9 +31,6 @@ python module f2hysop ! in\n
file(APPEND ${_file}
" ! Example
include '@CMAKE_SOURCE_DIR@/hysop/common_f/template.pyf'\n")
- file(APPEND ${_file}
- " ! precision
- !! include '@CMAKE_SOURCE_DIR@/hysop/f2py/parameters.pyf'\n")
if(WITH_FFTW)
file(APPEND ${_file}
" ! fftw
@@ -42,12 +39,12 @@ python module f2hysop ! in\n
if(WITH_SCALES)
file(APPEND ${_file}
" ! scales
- include '@CMAKE_SOURCE_DIR@/hysop/f2py/scales2py.pyf'\n")
+ include '@CMAKE_SOURCE_DIR@/hysop/scales_f/scales2py.pyf'\n")
endif()
if(WITH_EXTRAS)
file(APPEND ${_file}
" ! arnoldi
- include '@CMAKE_SOURCE_DIR@/hysop/fortran/arnoldi2py.pyf'\n")
+ include '@CMAKE_SOURCE_DIR@/hysop/numerics/extras_f/arnoldi.pyf'\n")
endif()
file(APPEND ${_file} " end interface\n
end python module f2hysop")
diff --git a/hysop/fortran/arnoldi2py.pyf b/hysop/fortran/arnoldi2py.pyf
deleted file mode 100644
index 6de3250e8665d2c260478dcd69baabd939ca45e4..0000000000000000000000000000000000000000
--- a/hysop/fortran/arnoldi2py.pyf
+++ /dev/null
@@ -1,16 +0,0 @@
-! -*- f90 -*-
-! Note: the context of this file is case sensitive.
-
-module arnoldi2py ! in arnoldi2py.f90
- use arnoldi
- use hysopparam
- use mpi
- subroutine arnoldi3d(solutions,ncli,ntot,nfp,nmodes,tps) ! in arnoldi2py.f90:arnoldi2py
- real(kind=pk) dimension(:,:),intent(inout) :: solutions
- integer intent(in) :: ncli
- integer intent(in) :: ntot
- integer intent(in) :: nfp
- integer intent(in) :: nmodes
- real(kind=pk) intent(in) :: tps
- end subroutine arnoldi3d
-end module arnoldi2py
diff --git a/hysop/fortran/template.f95 b/hysop/fortran/template.f95
deleted file mode 100644
index 1e5eadf578b0470214613a3e23dfaad8030d0aa3..0000000000000000000000000000000000000000
--- a/hysop/fortran/template.f95
+++ /dev/null
@@ -1,27 +0,0 @@
- !! Template file - Provide an example
- !! to write f2py interface between fortran and python
-
-
-module template_f2py
-
- implicit none
-
- !> a global var
- integer, parameter :: var1 = 12
-
-contains
-
- !> do something ...
- subroutine check_f2py(input, output)
- use precision
- !> input array
- real(kind=wp), dimension(:,:), intent(in) :: input
- !> output array
- real(kind=wp), dimension(:,:), intent(inout) :: output
-
- print *, 'template f2py for tab'
- output(:,:) = 2 * input(:,:)
- print *, 'aha hah', output(1,1)
- end subroutine check_f2py
-
-end module template_f2py
diff --git a/hysop/fortran/template.pyf b/hysop/fortran/template.pyf
deleted file mode 100644
index 301c7fa721a38a5834ad8c641d98240edb64838f..0000000000000000000000000000000000000000
--- a/hysop/fortran/template.pyf
+++ /dev/null
@@ -1,14 +0,0 @@
-! -*- f90 -*-
-! Note: the context of this file is case sensitive.
-
-module template_f2py ! in template.f95
- integer, parameter,optional :: var1=12
- subroutine check_f2py(input,output) ! in template.f95:template_f2py
- use precision
- real(kind=wp), dimension(:, :), intent(in) :: input
- real(kind=wp), dimension(size(input,1), size(input,2)), intent(in,out), depend(input,input) :: output
- end subroutine check_f2py
-end module template_f2py
-
-! This file was auto-generated with f2py (version:2).
-! See http://cens.ioc.ee/projects/f2py2e/
diff --git a/hysop/numerics/extras_f/arnoldi.f95 b/hysop/numerics/extras_f/arnoldi.f95
new file mode 100644
index 0000000000000000000000000000000000000000..0fef4f5fca58717d404424588688dff57cbf7d60
--- /dev/null
+++ b/hysop/numerics/extras_f/arnoldi.f95
@@ -0,0 +1,220 @@
+!=========================================================================
+!Computation of global modes using a time stepping (snapshots) technique
+!=========================================================================
+!=========================================================================
+!Reads snapshots, constructs the Hessenberg matrix
+!and computes the eigen-values/eigen-functions
+!=========================================================================
+module arnoldi
+
+ use precision
+ use parameters
+
+ implicit none
+
+contains
+
+ !======================
+ subroutine arnoldi3d(Mu,ncli,nt,nfp,nmodes,Tps)
+ implicit none
+ integer, intent(in) :: ncli ! number of snapshot
+ integer, intent(in) :: nt ! total number of points per snapshot
+ integer, intent(in) :: nmodes ! number of desired modes
+ integer, intent(in) :: nfp ! number of desired eigen functions
+ real(wp), intent(in) :: Tps ! sampling time step
+ real(wp), dimension(:,:), intent(inout) :: Mu ! snapshots
+
+ real(wp), dimension(:,:), allocatable :: un ! orthonomalized Krylov basis
+ real(wp), dimension(:,:), allocatable :: Hessenberg ! Hessenberg matrix
+ complex(wp), dimension(:), allocatable :: VP ! eigenvalues
+ complex(wp), dimension(:,:), allocatable :: FP, FP_J ! eigen functions
+ real(wp), dimension(:), allocatable :: reslog, res ! residuals
+ integer, dimension(:), allocatable :: t ! sorting array
+ real(wp), dimension(:), allocatable :: tab !
+ real(wp) :: norm, prod, error !
+
+ integer :: i,j,nmax
+
+ allocate(un(nt,ncli),Hessenberg(ncli,ncli-1),VP(ncli-1),FP(ncli-1,ncli-1),FP_J(nt,ncli))
+
+ nmax=0
+ VP=(0.,0.); FP=(0.,0.); FP_J=(0.,0.)
+ Hessenberg=0.
+
+ !==================================
+ ! Arnoldi method
+ !==================================
+
+ norm = dot_product(Mu(:,1),Mu(:,1))
+ norm = sqrt(norm)
+ un(:,1) = Mu(:,1)/norm ! first normalized vector u1
+
+ do j=2,ncli ! construct a normalized base u2... un
+ un(:,j)=Mu(:,j) ! w=Mu_j (We have Mu_j=U(:,j+1))
+ do i=1,j-1
+ Hessenberg(i,j-1)=dot_product(un(:,i),un(:,j))
+ un(:,j)=un(:,j)-un(:,i)*Hessenberg(i,j-1)
+ enddo
+
+ norm = dot_product(un(:,j),un(:,j))
+ Hessenberg(j,j-1) = sqrt(norm)
+
+ un(:,j) = un(:,j)/Hessenberg(j,j-1)! normalization
+
+ enddo
+
+ ! do i=1,nt
+ ! print *, 'Krylov basis:', un(i,:)
+ ! enddo
+
+ do i=1,ncli-1
+ print *, 'Hessenberg matrix:', Hessenberg(i,:)
+ enddo
+
+
+ !Check orthonormalization
+ !==================================
+
+ print *,'check ortho'
+
+ prod=0.
+ do i=1,ncli
+ do j=1,ncli
+ prod=dot_product(un(:,j),un(:,i))
+ if (abs(prod).gt.1e-14) then
+ print *,i,j,abs(prod)
+ endif
+ enddo
+ enddo
+
+
+ !Eigen-values and Eigen-functions related to Hessenberg matrix
+ ! +
+ !Eigen-values related to Jacobian matrix ==> spectra
+ !==============================================================
+
+ open(unit=10, file='spectrum.dat')
+ open(unit=11, file='spectrum_log.dat')
+
+ call spectrum(Hessenberg(1:ncli-1,:),ncli-1,VP,FP)
+
+ do i=1,ncli-1
+ write(10,*) dble(VP(i)), aimag(VP(i))
+ write(11,*) dble(log(VP(i)))/Tps, ATAN(aimag(VP(i)/DBLE(VP(i))))/Tps
+ enddo
+ close(10)
+ close(11)
+
+ !Eigen-functions related to Jacobian matrix
+ !==========================================
+ FP_J(1:nt,1:ncli-1)=matmul(un(1:nt,1:ncli-1),FP(1:ncli-1,1:ncli-1))
+ ! do i=1,nt
+ ! print *, 'FP_J', (FP_J(i,j),j=1,ncli-1)
+ ! enddo
+
+ !Residual calculation with respect to each mode
+ !==============================================
+
+ allocate(res(ncli-1),reslog(ncli-1))
+ error = Hessenberg(ncli,ncli-1)
+ print *,'last Hess',Hessenberg(ncli,ncli-1)
+
+ do i=1,ncli-1
+ res(i) = abs(FP(ncli-1,i))*error
+ reslog(i)=-log10(res(i))
+ print *,'residual',reslog(i),res(i)
+ enddo
+
+
+ !Modes are sorted with respect to residual
+ !==========================================
+ allocate(t(ncli-1))
+
+ do i=1,ncli-1
+ t(i)=i
+ enddo
+
+ call sort(reslog,ncli-1,t)
+
+ open(unit=201,file='spectrum_res.dat')
+ write(201,*)'VARIABLES ="WR","WI","RES"'
+
+ do i=1,nmodes
+ write(201,100) dble(log(VP(t(i))))/Tps,&
+ ATAN(aimag(VP(t(i))/DBLE(VP(t(i)))))/Tps,&
+ res(t(i))
+ enddo
+ close(201)
+ !
+ !Write the associated eigen functions
+ !====================================
+ ! allocate(tab(nfp))
+ !
+ ! open(unit=107, file='table.dat')
+ ! open(unit=108, file='spectrum_sorted.dat')
+ !
+ ! do i=1,nfp
+ !! call ecriture(FP_J(:,t(h)))
+ ! write(108,*) ATAN(aimag(VP(t(i))/DBLE(VP(t(i)))))/Tps,&
+ ! dble(log(VP(t(i))))/Tps
+ ! enddo
+ ! close(108)
+ !
+100 format (5(2x,e19.13))
+
+ end subroutine arnoldi3d
+
+ !===================
+ !Spectrum subroutine
+ !===================
+ subroutine spectrum(A,n,VP,VR)
+ implicit none
+ integer :: INFO
+ integer :: n,LWORK
+ real(wp), dimension(n,n) :: A
+ real(wp), dimension(:), allocatable :: RWORK
+ complex(wp), dimension(1,n) :: VL
+ complex(wp), dimension(n,n) :: VR
+ complex(wp), dimension(:), allocatable :: WORK
+ complex(wp), dimension(n):: VP
+
+ LWORK=4*n
+ allocate(WORK(LWORK),RWORK(2*n))
+ call ZGEEV('N','V', n, A*(1.,0.), n, VP, VL, 1, VR, n,&
+ WORK, LWORK, RWORK, INFO )
+ print *, 'VP', VP
+
+ end subroutine spectrum
+
+ !==================
+ !Sorting subroutine
+ !==================
+ subroutine sort(t,n,ind)
+ implicit none
+ integer :: i, j, n, tp1
+ real(wp), dimension(1:n) :: t
+ real(wp) :: temp
+ integer, dimension(1:n) :: ind
+
+ do i=1,n
+ do j=i+1,n
+ if ((t(i))<(t(j))) then
+
+ temp=t(i)
+ tp1=ind(i)
+
+ t(i)=t(j)
+ ind(i)=ind(j)
+
+ t(j)=temp
+ ind(j)=tp1
+
+ endif
+ enddo
+ enddo
+
+ return
+
+ end subroutine sort
+
+end module arnoldi
diff --git a/hysop/numerics/extras_f/arnoldi2py.pyf b/hysop/numerics/extras_f/arnoldi2py.pyf
new file mode 100644
index 0000000000000000000000000000000000000000..8dc28806a687c727c01df2d209f79776087fda23
--- /dev/null
+++ b/hysop/numerics/extras_f/arnoldi2py.pyf
@@ -0,0 +1,14 @@
+! -*- f90 -*-
+! Note: the context of this file is case sensitive.
+
+module arnoldi ! in arnoldi2py.f90
+ use precision
+ subroutine arnoldi3d(Mu,ncli,nt,nfp,nmodes,Tps) ! in arnoldi2py.f90:arnoldi2py
+ real(kind=wp) dimension(:,:),intent(inout) :: Mu
+ integer intent(in) :: ncli
+ integer intent(in) :: nt
+ integer intent(in) :: nfp
+ integer intent(in) :: nmodes
+ real(kind=pk) intent(in) :: Tps
+ end subroutine arnoldi3d
+end module arnoldi
diff --git a/hysop/operator/discrete/scales_advection.py b/hysop/operator/discrete/scales_advection.py
index 96ee6d4acb1980d44586513a76637c7a611246c2..8984877b805ddc7eb391be81ab4111e929391eb0 100644
--- a/hysop/operator/discrete/scales_advection.py
+++ b/hysop/operator/discrete/scales_advection.py
@@ -5,7 +5,9 @@
try:
from hysop.f2hysop import scales2py
except ImportError:
- from hysop.fakef2py import scales2py
+ msg = 'scales package not available for your hysop install.'
+ msg += 'Try to recompile with WITH_SCALES=ON'
+ raise ImportError(msg)
from hysop.operator.discrete.discrete import DiscreteOperator
from hysop.methods_keys import MultiScale
from hysop.constants import debug
diff --git a/hysop/operator/discrete/spectrum.py b/hysop/operator/discrete/spectrum.py
index 2bc03cbd5743b17753b18d8c1d3717523564baf7..a642770349c4b6774e578ed6a38fe7406624acf9 100755
--- a/hysop/operator/discrete/spectrum.py
+++ b/hysop/operator/discrete/spectrum.py
@@ -4,14 +4,15 @@
from hysop.operator.discrete.discrete import DiscreteOperator
from hysop.constants import debug, np, HYSOP_MPI_REAL
from hysop.tools.profiler import profile
-from hysop.tools.io_utils import IO
import hysop.tools.numpywrappers as npw
from hysop.mpi import MPI
-import os
-from hysop import __FFTW_ENABLED__
-if __FFTW_ENABLED__:
+try:
from hysop.f2hysop import fftw2py
+except ImportError:
+ msgE = 'fftw package not available for your hysop install.'
+ msgE += 'Try to recompile with WITH_FFTW=ON'
+ raise ImportError(msgE)
class FFTSpectrum(DiscreteOperator):
diff --git a/hysop/operator/tests/test_diff_poisson_3D.py b/hysop/operator/tests/test_diff_poisson_3D.py
index 0268a81acc2bac6be3205253dcccdc8364a51210..e406d2b4a6ef52f11a90bd2847047ee68e9fb1c2 100755
--- a/hysop/operator/tests/test_diff_poisson_3D.py
+++ b/hysop/operator/tests/test_diff_poisson_3D.py
@@ -1,10 +1,12 @@
-# -*- coding: utf-8 -*-
+"""Diffusion - Poisson sequence
+"""
import hysop as pp
from hysop.operator.poisson import Poisson
from hysop.operator.diffusion import Diffusion
from math import sqrt, pi, exp
from hysop.problem.simulation import Simulation
from hysop import testsenv
+from hysop.tools.parameters import Discretization
def computeVel(x, y, z):
@@ -38,12 +40,13 @@ def computeVort(x, y, z):
@testsenv.fftw_failed
-def test_Diff_Poisson():
+def test_diff_poisson():
+ """
+ """
# Parameters
nb = 33
boxLength = [1., 1., 1.]
boxMin = [0., 0., 0.]
- from hysop.tools.parameters import Discretization
d3D = Discretization([nb, nb, nb])
# Domain
@@ -56,7 +59,7 @@ def test_Diff_Poisson():
name='Vorticity', is_vector=True)
# FFT Diffusion operators and FFT Poisson solver
- diffusion = Diffusion(variables={vorti: d3D}, viscosity=0.002)
+ diffusion = Diffusion(vorticity=vorti, viscosity=0.002, discretization=d3D)
poisson = Poisson(velo, vorti, discretization=d3D)
diffusion.discretize()
diff --git a/hysop/operator/tests/test_spectrum.py b/hysop/operator/tests/test_spectrum.py
index e2404022e8f9f0f7117e1ef3c0cbdfa4492a63d4..cbd01d90c52484d05d4ea3d17413c12f34457fd2 100755
--- a/hysop/operator/tests/test_spectrum.py
+++ b/hysop/operator/tests/test_spectrum.py
@@ -35,5 +35,3 @@ def test_spectrum():
simu = Simulation(nb_iter=1)
simu.initialize()
op.apply(simu)
-
-
diff --git a/hysop/f2py/scales2py.f90 b/hysop/scales_f/scales2py.f90
similarity index 83%
rename from hysop/f2py/scales2py.f90
rename to hysop/scales_f/scales2py.f90
index c021618b45a1f723b4ea2176bf7e22036df5956b..165c071c14b68ab7be6a0e20576d73bd5d72524f 100755
--- a/hysop/f2py/scales2py.f90
+++ b/hysop/scales_f/scales2py.f90
@@ -6,8 +6,7 @@ use advec, only : advec_init,advec_step,advec_step_Inter_basic,advec_step_Inter_
use advec_vect, only : advec_step_Vect,advec_step_Inter_basic_Vect
use interpolation_velo, only : interpol_init
use mpi
-use hysopparam
-
+use precision
implicit none
@@ -23,14 +22,14 @@ contains
subroutine init_advection_solver(ncells,lengths,topodims,main_comm,datashape,offset,dim,order,dim_split)
integer, intent(in) :: dim
integer, dimension(dim),intent(in) :: ncells
- real(pk),dimension(dim), intent(in) :: lengths
+ real(wp),dimension(dim), intent(in) :: lengths
integer, dimension(dim), intent(in) :: topodims
integer, intent(in) :: main_comm
- integer(ik), dimension(dim), intent(out) :: datashape
- integer(ik), dimension(dim), intent(out) :: offset
+ integer(ip), dimension(dim), intent(out) :: datashape
+ integer(ip), dimension(dim), intent(out) :: offset
character(len=*), optional, intent(in) :: order, dim_split
- !! real(pk), optional, intent(out) :: stab_coeff
- real(pk) :: stab_coeff
+ !! real(wp), optional, intent(out) :: stab_coeff
+ real(wp) :: stab_coeff
!f2py integer optional , depend(ncells) :: dim=len(ncells)
!f2py intent(hide) dim
!f2py character(*) optional, intent(in) :: order = 'p_O2'
@@ -88,12 +87,12 @@ contains
!! @param[in,out] scal 3d scalar field which is advected
subroutine solve_advection(dt,vx,vy,vz,scal)
- real(pk), intent(in) :: dt
- real(pk), dimension(:,:,:), intent(in) :: vx, vy, vz
- real(pk), dimension(size(vx,1),size(vx,2),size(vx,3)), intent(inout) :: scal
- !f2py real(pk) intent(in,out), depend(size(vx,1)) :: scal
+ real(wp), intent(in) :: dt
+ real(wp), dimension(:,:,:), intent(in) :: vx, vy, vz
+ real(wp), dimension(size(vx,1),size(vx,2),size(vx,3)), intent(inout) :: scal
+ !f2py real(wp) intent(in,out), depend(size(vx,1)) :: scal
- real(pk) :: t0
+ real(wp) :: t0
t0 = MPI_Wtime()
call advec_step(dt,vx,vy,vz,scal)
@@ -110,12 +109,12 @@ contains
!! @param[in,out] cz 3d scalar field which is advected
subroutine solve_advection_vect(dt,vx,vy,vz,cx,cy,cz)
- real(pk), intent(in) :: dt
- real(pk), dimension(:,:,:), intent(in) :: vx, vy, vz
- real(pk), dimension(size(vx,1),size(vx,2),size(vx,3)), intent(inout) :: cx,cy,cz
- !f2py real(pk) intent(in,out), depend(size(vx,1)) :: cx
- !f2py real(pk) intent(in,out), depend(size(vx,1)) :: cy
- !f2py real(pk) intent(in,out), depend(size(vx,1)) :: cz
+ real(wp), intent(in) :: dt
+ real(wp), dimension(:,:,:), intent(in) :: vx, vy, vz
+ real(wp), dimension(size(vx,1),size(vx,2),size(vx,3)), intent(inout) :: cx,cy,cz
+ !f2py real(wp) intent(in,out), depend(size(vx,1)) :: cx
+ !f2py real(wp) intent(in,out), depend(size(vx,1)) :: cy
+ !f2py real(wp) intent(in,out), depend(size(vx,1)) :: cz
real(8) :: t0
@@ -134,9 +133,9 @@ contains
subroutine solve_advection_inter_basic(dt, Vx, Vy, Vz, scal)
! Input/Output
- real(pk), intent(in) :: dt
- real(pk), dimension(:,:,:), intent(in) :: Vx, Vy, Vz
- real(pk), dimension(:,:,:), intent(inout) :: scal
+ real(wp), intent(in) :: dt
+ real(wp), dimension(:,:,:), intent(in) :: Vx, Vy, Vz
+ real(wp), dimension(:,:,:), intent(inout) :: scal
!f2py intent(in,out) :: scal
call advec_step_Inter_basic(dt, Vx, Vy, Vz, scal)
@@ -154,9 +153,9 @@ contains
subroutine solve_advection_inter_basic_vec(dt, Vx, Vy, Vz, cx, cy, cz)
! Input/Output
- real(pk), intent(in) :: dt
- real(pk), dimension(:,:,:), intent(in) :: Vx, Vy, Vz
- real(pk), dimension(:,:,:), intent(inout) :: cx,cy,cz
+ real(wp), intent(in) :: dt
+ real(wp), dimension(:,:,:), intent(in) :: Vx, Vy, Vz
+ real(wp), dimension(:,:,:), intent(inout) :: cx,cy,cz
!f2py intent(in,out) :: cx
!f2py intent(in,out) :: cy
!f2py intent(in,out) :: cz
diff --git a/hysop/f2py/scales2py.pyf b/hysop/scales_f/scales2py.pyf
similarity index 57%
rename from hysop/f2py/scales2py.pyf
rename to hysop/scales_f/scales2py.pyf
index 584336cfac437a80667dccfa9dcc3cf15bcf2e87..d47582692be664ff51ae6838b391016281e914f5 100644
--- a/hysop/f2py/scales2py.pyf
+++ b/hysop/scales_f/scales2py.pyf
@@ -7,14 +7,14 @@ module scales2py ! in scales2py.f90
use cart_topology, only: cart_rank,coord,n_proc,discretisation_set_mesh_velo,cart_create,set_group_size,discretisation_create
use advec_vect, only: advec_step_vect,advec_step_inter_basic_vect
use mpi
- use hysopparam
+ use precision
subroutine init_advection_solver(ncells,lengths,topodims,main_comm,datashape,offset,dim,order,dim_split) ! in scales2py.f90:scales2py
integer dimension(dim),intent(in) :: ncells
- real(kind=pk) dimension(dim),intent(in),depend(dim) :: lengths
+ real(kind=wp) dimension(dim),intent(in),depend(dim) :: lengths
integer dimension(dim),intent(in),depend(dim) :: topodims
integer intent(in) :: main_comm
- integer(kind=ik) dimension(dim),intent(out),depend(dim) :: datashape
- integer(kind=ik) dimension(dim),intent(out),depend(dim) :: offset
+ integer(kind=ip) dimension(dim),intent(out),depend(dim) :: datashape
+ integer(kind=ip) dimension(dim),intent(out),depend(dim) :: offset
integer, optional,intent(hide),depend(ncells) :: dim=len(ncells)
character*(*), optional,intent(in) :: order='p_o2'
character*(*), optional,intent(in) :: dim_split
@@ -26,36 +26,36 @@ module scales2py ! in scales2py.f90
character*(*), optional,intent(in) :: formula
end subroutine init_multiscale
subroutine solve_advection(dt,vx,vy,vz,scal) ! in scales2py.f90:scales2py
- real(kind=pk) intent(in) :: dt
- real(kind=pk) dimension(:,:,:),intent(in) :: vx
- real(kind=pk) dimension(:,:,:),intent(in) :: vy
- real(kind=pk) dimension(:,:,:),intent(in) :: vz
- real(kind=pk) dimension(size(vx,1),size(vx,2),size(vx,3)),intent(in,out),depend(size(vx,1)) :: scal
+ real(kind=wp) intent(in) :: dt
+ real(kind=wp) dimension(:,:,:),intent(in) :: vx
+ real(kind=wp) dimension(:,:,:),intent(in) :: vy
+ real(kind=wp) dimension(:,:,:),intent(in) :: vz
+ real(kind=wp) dimension(size(vx,1),size(vx,2),size(vx,3)),intent(in,out),depend(size(vx,1)) :: scal
end subroutine solve_advection
subroutine solve_advection_vect(dt,vx,vy,vz,cx,cy,cz) ! in scales2py.f90:scales2py
- real(kind=pk) intent(in) :: dt
- real(kind=pk) dimension(:,:,:),intent(in) :: vx
- real(kind=pk) dimension(:,:,:),intent(in) :: vy
- real(kind=pk) dimension(:,:,:),intent(in) :: vz
- real(kind=pk) dimension(size(vx,1),size(vx,2),size(vx,3)),intent(in,out),depend(size(vx,1)) :: cx
- real(kind=pk) dimension(size(vx,1),size(vx,2),size(vx,3)),intent(in,out),depend(size(vx,1)) :: cy
- real(kind=pk) dimension(size(vx,1),size(vx,2),size(vx,3)),intent(in,out),depend(size(vx,1)) :: cz
+ real(kind=wp) intent(in) :: dt
+ real(kind=wp) dimension(:,:,:),intent(in) :: vx
+ real(kind=wp) dimension(:,:,:),intent(in) :: vy
+ real(kind=wp) dimension(:,:,:),intent(in) :: vz
+ real(kind=wp) dimension(size(vx,1),size(vx,2),size(vx,3)),intent(in,out),depend(size(vx,1)) :: cx
+ real(kind=wp) dimension(size(vx,1),size(vx,2),size(vx,3)),intent(in,out),depend(size(vx,1)) :: cy
+ real(kind=wp) dimension(size(vx,1),size(vx,2),size(vx,3)),intent(in,out),depend(size(vx,1)) :: cz
end subroutine solve_advection_vect
subroutine solve_advection_inter_basic(dt,vx,vy,vz,scal) ! in scales2py.f90:scales2py
- real(kind=pk) intent(in) :: dt
- real(kind=pk) dimension(:,:,:),intent(in) :: vx
- real(kind=pk) dimension(:,:,:),intent(in) :: vy
- real(kind=pk) dimension(:,:,:),intent(in) :: vz
- real(kind=pk) dimension(:,:,:),intent(in,out) :: scal
+ real(kind=wp) intent(in) :: dt
+ real(kind=wp) dimension(:,:,:),intent(in) :: vx
+ real(kind=wp) dimension(:,:,:),intent(in) :: vy
+ real(kind=wp) dimension(:,:,:),intent(in) :: vz
+ real(kind=wp) dimension(:,:,:),intent(in,out) :: scal
end subroutine solve_advection_inter_basic
subroutine solve_advection_inter_basic_vec(dt,vx,vy,vz,cx,cy,cz) ! in scales2py.f90:scales2py
- real(kind=pk) intent(in) :: dt
- real(kind=pk) dimension(:,:,:),intent(in) :: vx
- real(kind=pk) dimension(:,:,:),intent(in) :: vy
- real(kind=pk) dimension(:,:,:),intent(in) :: vz
- real(kind=pk) dimension(:,:,:),intent(in,out) :: cx
- real(kind=pk) dimension(:,:,:),intent(in,out) :: cy
- real(kind=pk) dimension(:,:,:),intent(in,out) :: cz
+ real(kind=wp) intent(in) :: dt
+ real(kind=wp) dimension(:,:,:),intent(in) :: vx
+ real(kind=wp) dimension(:,:,:),intent(in) :: vy
+ real(kind=wp) dimension(:,:,:),intent(in) :: vz
+ real(kind=wp) dimension(:,:,:),intent(in,out) :: cx
+ real(kind=wp) dimension(:,:,:),intent(in,out) :: cy
+ real(kind=wp) dimension(:,:,:),intent(in,out) :: cz
end subroutine solve_advection_inter_basic_vec
end module scales2py
diff --git a/setup.py.in b/setup.py.in
index f7c15653caa3b1cafec0b71c281063dba73e77f1..9fa5158b74f458ba6a9d7d4991efd0955472aee4 100644
--- a/setup.py.in
+++ b/setup.py.in
@@ -70,6 +70,18 @@ def parseCMakeDefines(var):
# --- external libraries to be linked with ---
hysop_link_libraries = parseCMakeVar("@HYSOP_LINK_LIBRARIES@")
+
+# deal with macosx framework ...
+extra_flags = []
+for lib in hysop_link_libraries:
+ res = lib.find('framework')
+ if res >= 0:
+ index = hysop_link_libraries.index(lib)
+ filename = os.path.basename(lib)
+ hysop_link_libraries.pop(index)
+ filename = filename.split('.')[0]
+ extra_flags.append('-framework ' + filename)
+hysop_link_libraries += extra_flags
#hysop_link_libraries_names = set([])
#hysop_link_libraries_dirs = set([])
# use set to avoid dupl.
@@ -266,8 +278,8 @@ ext = {}
if enable_fortran is "ON":
fortran_root = \
'@CMAKE_SOURCE_DIR@/hysop'
- hysop_libdir = [] ##### ['@CMAKE_BINARY_DIR@/src']
- #####hysoplib = ['@HYSOP_LIBRARY_NAME@']
+ hysop_libdir = ['@CMAKE_BINARY_DIR@/src']
+ hysoplib = ['@HYSOP_LIBRARY_NAME@']
f2py_options = ['--no-lower']
fortran_src = set([])
# -- list of directories which contains fortran sources --
@@ -278,18 +290,21 @@ if enable_fortran is "ON":
# -- fftw fortran sources --
withfftw = "@WITH_FFTW@"
if withfftw is "ON":
- # fortran_src.add('f2py/fftw2py.f90')
- subdirs.append(os.path.join('numerics/fftw_f'))
- #fortran_src.add('f2py/fftw2py.f90')
+ subdirs.append(os.path.join('numerics', 'fftw_f'))
fftwdir = '@FFTWLIB@'
- #hysoplib.append('fftw3')
- #hysoplib.append('fftw3_mpi')
hysop_libdir.append(fftwdir)
# -- scales sources --
withscales = '@WITH_SCALES@'
if withscales is "ON":
- fortran_src.add('f2py/scales2py.f90')
+ subdirs.append('scales_f')
+ #fortran_src.add('f2py/scales2py.f90')
+
+
+ # -- other fortran sources --
+ withextras = '@WITH_EXTRAS@'
+ if withextras is "ON":
+ subdirs.append(os.path.join('numerics', 'extras_f'))
# -- set full path to fortran sources --
fortran_src = list(fortran_src)
@@ -316,8 +331,8 @@ if enable_fortran is "ON":
ext['f2hysop'] = create_fortran_extension(
name='hysop.f2hysop',
sources=fortran_src,
- #libdir=hysop_libdir,
- #libs=hysoplib,
+ libdir=hysop_libdir,
+ libs=hysoplib,
pyf_file=pyf_file,
src_dirs=num_dirs)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a66936986a42225f07626ad70dc525877387296c..a3d17a69c40c655407f1087e02f7df41d14dd11a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -10,11 +10,6 @@
set(${HYSOP_LIBRARY_NAME}_SRCDIRS
.)
-# --- fftw interface ---
-if(WITH_FFTW)
- list(APPEND ${HYSOP_LIBRARY_NAME}_SRCDIRS fftw)
-endif()
-
if(USE_CXX)
list(APPEND ${HYSOP_LIBRARY_NAME}_SRCDIRS hysop++/src)
endif()
@@ -81,25 +76,25 @@ target_link_libraries(${HYSOP_LIBRARY_NAME} PRIVATE ${HYSOP_LINK_LIBRARIES})
# way that does not depends on python.
# At the time it only includes fftw tests.
if(WITH_FORTRAN_TESTS)
- if(WITH_FFTW)
- # Set the name of a executable file that will be linked with libHYSOP_LIBRARY_NAME.
- # Useful to test libhysop in a way that does not depend on python.
- set(EXE_NAME ${HYSOP_LIBRARY_NAME}Run)
- # A main file to create an executable (test purpose)
- # Any files in these dirs will be used to create HySoP exec (linked with libhysop)
- set(${EXE_NAME}_SRCDIRS main)
- # Source and header files list, to generate a working executable based on libhysop.
- get_sources("${${EXE_NAME}_SRCDIRS}")
- get_headers("${${EXE_NAME}_SRCDIRS}")
- set(${EXE_NAME}_SRCS ${SOURCES_FILES})
- set(${EXE_NAME}_HDRS ${HDRS_FILES})
- list(APPEND ${EXE_NAME}_SRC ${${EXE_NAME}_HDRS})
- include_directories(${${EXE_NAME}_HDRS})
- add_executable(${EXE_NAME} ${${EXE_NAME}_SRCS})
- add_dependencies(${EXE_NAME} ${HYSOP_LIBRARY_NAME})
+ # if(WITH_FFTW)
+ # # Set the name of a executable file that will be linked with libHYSOP_LIBRARY_NAME.
+ # # Useful to test libhysop in a way that does not depend on python.
+ # set(EXE_NAME ${HYSOP_LIBRARY_NAME}Run)
+ # # A main file to create an executable (test purpose)
+ # # Any files in these dirs will be used to create HySoP exec (linked with libhysop)
+ # set(${EXE_NAME}_SRCDIRS main)
+ # # Source and header files list, to generate a working executable based on libhysop.
+ # get_sources("${${EXE_NAME}_SRCDIRS}")
+ # get_headers("${${EXE_NAME}_SRCDIRS}")
+ # set(${EXE_NAME}_SRCS ${SOURCES_FILES})
+ # set(${EXE_NAME}_HDRS ${HDRS_FILES})
+ # list(APPEND ${EXE_NAME}_SRC ${${EXE_NAME}_HDRS})
+ # include_directories(${${EXE_NAME}_HDRS})
+ # add_executable(${EXE_NAME} ${${EXE_NAME}_SRCS})
+ # add_dependencies(${EXE_NAME} ${HYSOP_LIBRARY_NAME})
- # libs to link with EXE_NAME
- target_link_libraries(${EXE_NAME} ${HYSOP_LIBRARY_NAME})
- target_link_libraries(${EXE_NAME} ${HYSOP_LINK_LIBRARIES})
- endif()
+ # # libs to link with EXE_NAME
+ # target_link_libraries(${EXE_NAME} ${HYSOP_LIBRARY_NAME})
+ # target_link_libraries(${EXE_NAME} ${HYSOP_LINK_LIBRARIES})
+ # endif()
endif()
diff --git a/src/extras/arnoldi.f90 b/src/extras/arnoldi.f90
deleted file mode 100644
index 386498094fa218387456692073ffcafff4cf3ef2..0000000000000000000000000000000000000000
--- a/src/extras/arnoldi.f90
+++ /dev/null
@@ -1,219 +0,0 @@
-!=========================================================================
-!Computation of global modes using a time stepping (snapshots) technique
-!=========================================================================
-!=========================================================================
-!Reads snapshots, constructs the Hessenberg matrix
-!and computes the eigen-values/eigen-functions
-!=========================================================================
-module arnoldi
-
- use client_data
-
- implicit none
-
-contains
-
-!======================
-subroutine arnoldi3d(Mu,ncli,nt,nfp,nmodes,Tps)
- implicit none
- integer, intent(in) :: ncli ! number of snapshot
- integer, intent(in) :: nt ! total number of points per snapshot
- integer, intent(in) :: nmodes ! number of desired modes
- integer, intent(in) :: nfp ! number of desired eigen functions
- real(mk), intent(in) :: Tps ! sampling time step
- real(mk), dimension(:,:), intent(inout) :: Mu ! snapshots
-
- real(mk), dimension(:,:), allocatable :: un ! orthonomalized Krylov basis
- real(mk), dimension(:,:), allocatable :: Hessenberg ! Hessenberg matrix
- complex(mk), dimension(:), allocatable :: VP ! eigenvalues
- complex(mk), dimension(:,:), allocatable :: FP, FP_J ! eigen functions
- real(mk), dimension(:), allocatable :: reslog, res ! residuals
- integer, dimension(:), allocatable :: t ! sorting array
- real(mk), dimension(:), allocatable :: tab !
- real(mk) :: norm, prod, error !
-
- integer :: i,j,k,nmax
-
- allocate(un(nt,ncli),Hessenberg(ncli,ncli-1),VP(ncli-1),FP(ncli-1,ncli-1),FP_J(nt,ncli))
-
- nmax=0
- VP=(0.,0.); FP=(0.,0.); FP_J=(0.,0.)
- Hessenberg=0.
-
- !==================================
- ! Arnoldi method
- !==================================
-
- norm = dot_product(Mu(:,1),Mu(:,1))
- norm = sqrt(norm)
- un(:,1) = Mu(:,1)/norm ! first normalized vector u1
-
- do j=2,ncli ! construct a normalized base u2... un
- un(:,j)=Mu(:,j) ! w=Mu_j (We have Mu_j=U(:,j+1))
- do i=1,j-1
- Hessenberg(i,j-1)=dot_product(un(:,i),un(:,j))
- un(:,j)=un(:,j)-un(:,i)*Hessenberg(i,j-1)
- enddo
-
- norm = dot_product(un(:,j),un(:,j))
- Hessenberg(j,j-1) = sqrt(norm)
-
- un(:,j) = un(:,j)/Hessenberg(j,j-1)! normalization
-
- enddo
-
-! do i=1,nt
-! print *, 'Krylov basis:', un(i,:)
-! enddo
-
- do i=1,ncli-1
- print *, 'Hessenberg matrix:', Hessenberg(i,:)
- enddo
-
-
-!Check orthonormalization
-!==================================
-
- print *,'check ortho'
-
- prod=0.
- do i=1,ncli
- do j=1,ncli
- prod=dot_product(un(:,j),un(:,i))
- if (abs(prod).gt.1e-14) then
- print *,i,j,abs(prod)
- endif
- enddo
- enddo
-
-
-!Eigen-values and Eigen-functions related to Hessenberg matrix
-! +
-!Eigen-values related to Jacobian matrix ==> spectra
-!==============================================================
-
- open(unit=10, file='spectrum.dat')
- open(unit=11, file='spectrum_log.dat')
-
-call spectrum(Hessenberg(1:ncli-1,:),ncli-1,VP,FP)
-
- do i=1,ncli-1
- write(10,*) dble(VP(i)), aimag(VP(i))
- write(11,*) dble(log(VP(i)))/Tps, ATAN(aimag(VP(i)/DBLE(VP(i))))/Tps
- enddo
- close(10)
- close(11)
-
-!Eigen-functions related to Jacobian matrix
-!==========================================
- FP_J(1:nt,1:ncli-1)=matmul(un(1:nt,1:ncli-1),FP(1:ncli-1,1:ncli-1))
-! do i=1,nt
-! print *, 'FP_J', (FP_J(i,j),j=1,ncli-1)
-! enddo
-
-!Residual calculation with respect to each mode
-!==============================================
-
- allocate(res(ncli-1),reslog(ncli-1))
- error = Hessenberg(ncli,ncli-1)
- print *,'last Hess',Hessenberg(ncli,ncli-1)
-
- do i=1,ncli-1
- res(i) = abs(FP(ncli-1,i))*error
- reslog(i)=-log10(res(i))
- print *,'residual',reslog(i),res(i)
- enddo
-
-
-!Modes are sorted with respect to residual
-!==========================================
- allocate(t(ncli-1))
-
- do i=1,ncli-1
- t(i)=i
- enddo
-
- call sort(reslog,ncli-1,t)
-
- open(unit=201,file='spectrum_res.dat')
- write(201,*)'VARIABLES ="WR","WI","RES"'
-
- do i=1,nmodes
- write(201,100) dble(log(VP(t(i))))/Tps,&
- ATAN(aimag(VP(t(i))/DBLE(VP(t(i)))))/Tps,&
- res(t(i))
- enddo
- close(201)
-!
-!Write the associated eigen functions
-!====================================
-! allocate(tab(nfp))
-!
-! open(unit=107, file='table.dat')
-! open(unit=108, file='spectrum_sorted.dat')
-!
-! do i=1,nfp
-!! call ecriture(FP_J(:,t(h)))
-! write(108,*) ATAN(aimag(VP(t(i))/DBLE(VP(t(i)))))/Tps,&
-! dble(log(VP(t(i))))/Tps
-! enddo
-! close(108)
-!
-100 format (5(2x,e19.13))
-
-end subroutine arnoldi3d
-
-!===================
-!Spectrum subroutine
-!===================
-subroutine spectrum(A,n,VP,VR)
- implicit none
- integer :: INFO
- integer :: n,LWORK
- real(mk), dimension(n,n) :: A
- real(mk), dimension(:), allocatable :: RWORK
- complex(mk), dimension(1,n) :: VL
- complex(mk), dimension(n,n) :: VR
- complex(mk), dimension(:), allocatable :: WORK
- complex(mk), dimension(n):: VP
-
- LWORK=4*n
- allocate(WORK(LWORK),RWORK(2*n))
- call ZGEEV('N','V', n, A*(1.,0.), n, VP, VL, 1, VR, n,&
- WORK, LWORK, RWORK, INFO )
- print *, 'VP', VP
-
-end subroutine spectrum
-
-!==================
-!Sorting subroutine
-!==================
-subroutine sort(t,n,ind)
- implicit none
- integer :: i, j, n, tp1
- real(mk), dimension(1:n) :: t
- real(mk) :: temp
- integer, dimension(1:n) :: ind
-
- do i=1,n
- do j=i+1,n
- if ((t(i))<(t(j))) then
-
- temp=t(i)
- tp1=ind(i)
-
- t(i)=t(j)
- ind(i)=ind(j)
-
- t(j)=temp
- ind(j)=tp1
-
- endif
- enddo
- enddo
-
- return
-
-end subroutine sort
-
-end module arnoldi