Merge pull request #27 from mwiesenberger/x-point

X point geometries

README.adoc

@@ -307,28 +307,38 @@ version including informative pdf writeups on implemented numerical
 methods directly from source code. This depends on the `doxygen`,
 `libjs-mathjax` and `graphviz` packages and LateX. Type `make doc` in
 the folder `path/to/feltor/doc` and open `index.html` (a symbolic link
-to `dg/html/modules.html`) with your favorite browser. 
+to `dg/html/modules.html`) with your favorite browser.
 Finally, also note the documentations of https://thrust.github.io/doc/modules.html[thrust]
 and https://cusplibrary.github.io/[cusp].
 
 == 3. Acknowledgements
 
-FELTOR is developed by Matthias Wiesenberger and Markus Held and
-receives contributions from an increasing number of people. We
-gratefully acknowledge fruitful discussions and code contribution from
+FELTOR's main authors are Matthias Wiesenberger and Markus Held.
+We received funding from the Austrian Science Fund (FWF) within
+project W1227 and Y398 under project leader Alexander Kendl at
+Innsbruck University.
+The work has been carried out within
+the framework of the EUROfusion Consortium and has received funding
+from the Euratom research and training programme 2014-2018 under
+grant agreement No 633053.
+MW has received funding from the European Union’s Horizon 2020
+research and innovation programme under the Marie
+Sklodowska-Curie grant agreement no. 713683 (COFUNDfellowsDTU).
+
+We gratefully acknowledge fruitful discussions and code contribution from
 
 - Ralph Kube
-- Eduard Reiter 
-- Lukas Einkemmer 
+- Eduard Reiter
+- Lukas Einkemmer
 - Jakob Gath
 
 We further acknowledge support on the Knights landing architecture from
-the High Level Support Team from 
+the High Level Support Team from
 
-- Albert Gutiérrez 
+- Albert Gutiérrez
 - Xavier Saez
 
-and from Intel Barcelona 
+and from Intel Barcelona
 
 - Harald Servat
 

config/default.mk

@@ -4,10 +4,10 @@ INCLUDED=1
 #compiler and compiler options
 CC=g++ #C++ compiler
 MPICC=mpic++  #mpi compiler
-CFLAGS=-Wall -std=c++11  #flags for CC
+CFLAGS=-Wall -std=c++11 -mavx #flags for CC
 NVCC=nvcc #CUDA compiler
-NVCCARCH=-arch sm_20 #nvcc gpu compute capability
-NVCCFLAGS= -std=c++11 -Xcompiler -Wall#flags for NVCC
+NVCCARCH=-arch sm_35 -Xcudafe "--diag_suppress=code_is_unreachable --diag_suppress=initialization_not_reachable" #nvcc gpu compute capability
+NVCCFLAGS= -std=c++11 -Xcompiler "-Wall -mavx" # -mfma"#flags for NVCC
 OPT=-O3 # optimization flags for host code
 OMPFLAG=-fopenmp #openmp flag for CC and MPICC
 

config/marconi.mk

@@ -5,12 +5,12 @@ MPICC=mpiicc
 OPT=-O3 -xHost  # overwritten for mic in devices.mk
 #MPICFLAGS+= -DMPICH_IGNORE_CXX_SEEK
 OMPFLAG=-qopenmp
-CFLAGS+= -restrict
+CFLAGS=-Wall -std=c++11 -restrict #-mfma  #flags for CC
 
 INCLUDE += -I$(HOME)/include # cusp, thrust
 INCLUDE += -I$(NETCDF_INC) -I$(HDF5_INC)
 LIBS    +=-L$(HDF5_LIB) -lhdf5 -lhdf5_hl
-LIBS    +=-L$(NETCDF_LIB) -lnetcdf -lcurl
+LIBS    +=-L$(NETCDF_LIB) -lnetcdf #-lcurl
 GLFLAGS  = -lm 
 endif
 #############################modules to load in .bashrc#######################

diag/feltordiag.cu

@@ -43,7 +43,7 @@ int main( int argc, char* argv[])
     Json::Reader reader;
     Json::Value js,gs;
     reader.parse( input, js, false);
-    const eule::Parameters p(js);
+    const feltor::Parameters p(js);
     reader.parse( geom, gs, false);
     const dg::geo::solovev::Parameters gp(gs);
     p.display();

diag/filamentdiag.cu

@@ -38,7 +38,7 @@ int main( int argc, char* argv[])
     Json::Reader reader;
     Json::Value js,gs;
     reader.parse( input, js, false);
-    const eule::Parameters p(js);
+    const feltor::Parameters p(js);
     reader.parse( geom, gs, false);
     const dg::geo::solovev::Parameters gp(gs);
     p.display();

diag/ncdiag.cpp

@@ -44,7 +44,7 @@ int main( int argc, char* argv[])
     Json::Reader reader;
     Json::Value js,gs;
     reader.parse( input, js, false);
-    const eule::Parameters p(js);
+    const feltor::Parameters p(js);
     reader.parse( geom, gs, false);
     const dg::geo::solovev::Parameters gp(gs);
     p.display();

doc/header.html

@@ -25,7 +25,7 @@
   <ul class="tablist"> 
       <!-- <li><a href="../../../index.html">/</a><li>-->
     <li><a href="../../dg/html/modules.html">dG</a></li>
-    <li><a href="../../geometries/html/modules.html">geometry</a></li>
+    <li><a href="../../geometries/html/modules.html">geometries</a></li>
     <li><a href="../../file/html/namespacefile.html">file</a></li>
   </ul>
   <!--End My own title area-->

inc/dg/Doxyfile

@@ -1523,7 +1523,7 @@ MATHJAX_FORMAT         = HTML-CSS
 # The default value is: http://cdn.mathjax.org/mathjax/latest.
 # This tag requires that the tag USE_MATHJAX is set to YES.
 
-MATHJAX_RELPATH        = https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?... 
+MATHJAX_RELPATH        = https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?...
 
 # The MATHJAX_EXTENSIONS tag can be used to specify one or more MathJax
 # extension names that should be enabled during MathJax rendering. For example
@@ -2110,7 +2110,7 @@ TAGFILES               =
 # tag file that is based on the input files it reads. See section "Linking to
 # external documentation" for more information about the usage of tag files.
 
-GENERATE_TAGFILE       = 
+GENERATE_TAGFILE       =
 
 # If the ALLEXTERNALS tag is set to YES, all external class will be listed in
 # the class index. If set to NO, only the inherited external classes will be

inc/dg/Makefile

@@ -1,8 +1,8 @@
 device=gpu
 
-#configure machine 
+#configure machine
 include ../../config/default.mk
-include ../../config/*.mk 
+include ../../config/*.mk
 include ../../config/devices/devices.mk
 
 INCLUDE+= -I../../ # other project libraries
@@ -16,22 +16,22 @@ all: $(CPPFILES:%.cpp=%) $(CUFILES:%.cu=%)
 %_t: %_t.cu
 	$(CC) $(OPT) $(INCLUDE) -DDG_DEBUG $(CFLAGS) $< -o $@  -g
 
-%_b: %_b.cu 
-	$(CC) $(OPT) $(CFLAGS) -DDG_BENCHMARK $< -o $@ $(INCLUDE) -g 
+%_b: %_b.cu
+	$(CC) $(OPT) $(CFLAGS) -DDG_BENCHMARK $< -o $@ $(INCLUDE) -g
 
-%_mpit: %_mpit.cu 
+%_mpit: %_mpit.cu
 	$(MPICC) $(INCLUDE) -DDG_DEBUG $(MPICFLAGS) $< -o $@ -g
 
 %_mpib: %_mpib.cu
 	$(MPICC) $(OPT) $(MPICFLAGS) -DDG_BENCHMARK $< -o $@ $(INCLUDE) -g
 
-bathRZ_t: bathRZ_t.cu 
+bathRZ_t: bathRZ_t.cu
 	$(CC) $(OPT) $(CFLAGS) $< -o $@ $(GLFLAGS) $(INCLUDE)  -g
-	
+
 .PHONY: clean doc
 
-doc: 
+doc:
 	doxygen Doxyfile
 
 clean:
-	rm -f *_t *_b *_mpit *_mpib 
+	rm -f *_t *_b *_mpit *_mpib

inc/dg/README

@@ -1,2 +1,2 @@
-contains container free numerical algorithms based on the template metaprogramming approach and test programs 
+contains container free numerical algorithms based on the template metaprogramming approach and test programs
 

inc/dg/arakawa.h

@@ -11,8 +11,8 @@
 #include "backend/mpi_evaluation.h"
 #endif
 
-/*! @file 
-  
+/*! @file
+
   @brief object for computation of Poisson bracket
   */
 namespace dg
@@ -22,7 +22,7 @@ namespace dg
  * @brief X-space generalized version of Arakawa's scheme
  *
  * Computes \f[ [f,g] := 1/\sqrt{g_{2d}}\left(\partial_x f\partial_y g - \partial_y f\partial_x g\right) \f]
- * where \f$ g_{2d} = g/g_{zz}\f$ is the two-dimensional volume element of the plane in 2x1 product space. 
+ * where \f$ g_{2d} = g/g_{zz}\f$ is the two-dimensional volume element of the plane in 2x1 product space.
  * @snippet arakawa_t.cu function
  * @snippet arakawa_t.cu doxygen
  * @copydoc hide_geometry_matrix_container
@@ -48,15 +48,15 @@ struct ArakawaX
      * @brief Compute poisson's bracket (25+2 memops)
      *
      * Computes \f[ [f,g] := 1/\sqrt{g_{2d}}\left(\partial_x f\partial_y g - \partial_y f\partial_x g\right) \f]
-     * where \f$ g_{2d} = g/g_{zz}\f$ is the two-dimensional volume element of the plane in 2x1 product space. 
+     * where \f$ g_{2d} = g/g_{zz}\f$ is the two-dimensional volume element of the plane in 2x1 product space.
      * @param lhs left hand side in x-space
      * @param rhs rights hand side in x-space
      * @param result Poisson's bracket in x-space
      */
     void operator()( const container& lhs, const container& rhs, container& result);
 
     /**
-     * @brief Return internally used x - derivative 
+     * @brief Return internally used x - derivative
      *
      * The same as a call to dg::create::dx( g, bcx)
      * @return derivative
@@ -71,11 +71,11 @@ struct ArakawaX
     const Matrix& dy() {return bdyf;}
 
     /**
-     * @brief Compute the total variation integrand 
+     * @brief Compute the total variation integrand
      *
      * Computes \f[ (\nabla\phi)^2 = \partial_i \phi g^{ij}\partial_j \phi \f]
      * in the plane of a 2x1 product space
-     * @param phi function 
+     * @param phi function
      * @param varphi may equal phi, contains result on output
      */
     void variation( const container& phi, container& varphi)
@@ -95,8 +95,8 @@ struct ArakawaX
 };
 ///@cond
 template<class Geometry, class Matrix, class container>
-ArakawaX<Geometry, Matrix, container>::ArakawaX( const Geometry& g ): 
-    dxlhs( dg::evaluate( one, g) ), dxrhs(dxlhs), dylhs(dxlhs), dyrhs( dxlhs), helper_( dxlhs), 
+ArakawaX<Geometry, Matrix, container>::ArakawaX( const Geometry& g ):
+    dxlhs( dg::evaluate( one, g) ), dxrhs(dxlhs), dylhs(dxlhs), dyrhs( dxlhs), helper_( dxlhs),
     bdxf( dg::create::dx( g, g.bcx())),
     bdyf( dg::create::dy( g, g.bcy()))
 {
@@ -105,11 +105,11 @@ ArakawaX<Geometry, Matrix, container>::ArakawaX( const Geometry& g ):
     dg::tensor::sqrt(perp_vol_inv_);
 }
 template<class Geometry, class Matrix, class container>
-ArakawaX<Geometry, Matrix, container>::ArakawaX( const Geometry& g, bc bcx, bc bcy): 
+ArakawaX<Geometry, Matrix, container>::ArakawaX( const Geometry& g, bc bcx, bc bcy):
     dxlhs( dg::evaluate( one, g) ), dxrhs(dxlhs), dylhs(dxlhs), dyrhs( dxlhs), helper_( dxlhs),
     bdxf(dg::create::dx( g, bcx)),
     bdyf(dg::create::dy( g, bcy))
-{ 
+{
     metric_=g.metric().perp();
     perp_vol_inv_ = dg::tensor::determinant(metric_);
     dg::tensor::sqrt(perp_vol_inv_);

inc/dg/arakawa3d_b.cu

@@ -9,18 +9,18 @@
 const double lx = 2.*M_PI;
 const double ly = 2.*M_PI;
 double left( double x, double y) {return sin(x)*sin(y);}
-double right( double x, double y) {return sin(2*x)*sin(2*y);} 
+double right( double x, double y) {return sin(2*x)*sin(2*y);}
 //double right2( double x, double y) {return sin(y);}
-double jacobian( double x, double y) 
+double jacobian( double x, double y)
 {
     return cos(x)*sin(y)*2*sin(2*x)*cos(2*y)-sin(x)*cos(y)*2*cos(2*x)*sin(2*y);
 }
 
 const double lz = 1.;
 double left( double x, double y, double z) {return left(x,y)*z;}
-double right( double x, double y, double z) {return right(x,y)*z;} 
+double right( double x, double y, double z) {return right(x,y)*z;}
 //double right2( double x, double y) {return sin(y);}
-double jacobian( double x, double y, double z) 
+double jacobian( double x, double y, double z)
 {
     return jacobian(x,y)*z*z;
 }
@@ -44,7 +44,7 @@ int main()
     dg::Timer t;
 
     dg::ArakawaX<dg::CartesianGrid3d, dg::DMatrix, dg::DVec> arakawa( grid);
-    t.tic(); 
+    t.tic();
     for( unsigned i=0; i<20; i++)
         arakawa( lhs, rhs, jac);
     t.toc();
@@ -56,7 +56,7 @@ int main()
     std::cout << "Mean rhs*Jacobian is "<<dg::blas2::dot( rhs, w3d, jac)<<"\n";
     std::cout << "Mean   n*Jacobian is "<<dg::blas2::dot( lhs, w3d, jac)<<"\n";
     dg::blas1::axpby( 1., sol, -1., jac);
-    std::cout << "Distance to solution "<<sqrt(dg::blas2::dot( w3d, jac))<<std::endl; 
+    std::cout << "Distance to solution "<<sqrt(dg::blas2::dot( w3d, jac))<<std::endl;
     return 0;
 }
 

inc/dg/arakawa3d_mpib.cu

@@ -20,25 +20,25 @@ const double ly = 2*M_PI;
 /*
 double left( double x, double y) { return sin(x)*cos(y);}
 double right( double x, double y){ return exp(0.1*(x+y)); }
-double jacobian( double x, double y) 
+double jacobian( double x, double y)
 {
-    return exp( x-M_PI)*(sin(x)+cos(x))*sin(y) * exp(y-M_PI)*sin(x)*(sin(y) + cos(y)) - sin(x)*exp(x-M_PI)*cos(y) * cos(x)*sin(y)*exp(y-M_PI); 
+    return exp( x-M_PI)*(sin(x)+cos(x))*sin(y) * exp(y-M_PI)*sin(x)*(sin(y) + cos(y)) - sin(x)*exp(x-M_PI)*cos(y) * cos(x)*sin(y)*exp(y-M_PI);
 }
 */
 
 dg::bc bcx = dg::PER;
 dg::bc bcy = dg::PER;
 double left( double x, double y) {return sin(x)*cos(y);}
 double right( double x, double y) {return cos(x)*sin(y);}
-double jacobian( double x, double y) 
+double jacobian( double x, double y)
 {
-    return cos(x)*cos(y)*cos(x)*cos(y) - sin(x)*sin(y)*sin(x)*sin(y); 
+    return cos(x)*cos(y)*cos(x)*cos(y) - sin(x)*sin(y)*sin(x)*sin(y);
 }
 const double lz = 1.;
 double left( double x, double y, double z) {return left(x,y)*z;}
-double right( double x, double y, double z) {return right(x,y)*z;} 
+double right( double x, double y, double z) {return right(x,y)*z;}
 //double right2( double x, double y) {return sin(y);}
-double jacobian( double x, double y, double z) 
+double jacobian( double x, double y, double z)
 {
     return jacobian(x,y)*z*z;
 }
@@ -49,7 +49,7 @@ int main(int argc, char* argv[])
 {
     MPI_Init( &argc, &argv);
     int rank;
-    unsigned n, Nx, Ny, Nz; 
+    unsigned n, Nx, Ny, Nz;
     MPI_Comm comm;
     mpi_init3d( bcx, bcy,dg::PER, n, Nx, Ny,Nz, comm);
     dg::MPIGrid3d grid( 0, lx, 0, ly, 0,lz, n, Nx, Ny, Nz, bcx, bcy, dg::PER, comm);
@@ -64,7 +64,7 @@ int main(int argc, char* argv[])
 
     dg::ArakawaX<dg::CartesianMPIGrid3d, Matrix, Vector> arakawa( grid);
     unsigned multi=20;
-    t.tic(); 
+    t.tic();
     for( unsigned i=0; i<multi; i++)
         arakawa( lhs, rhs, jac);
     t.toc();
@@ -80,7 +80,7 @@ int main(int argc, char* argv[])
     if(rank==0) std::cout << "Mean lhs*Jacobian is "<<result<<"\n";
     dg::blas1::axpby( 1., sol, -1., jac);
     result = sqrt( dg::blas2::dot( w3d, jac));
-    if(rank==0) std::cout << "Distance to solution "<<result<<std::endl; 
+    if(rank==0) std::cout << "Distance to solution "<<result<<std::endl;
 
 
     MPI_Finalize();

inc/dg/arakawa_b.cu

@@ -21,19 +21,19 @@ const double ly = 2*M_PI;
 /*
 double left( double x, double y) { return sin(x)*cos(y);}
 double right( double x, double y){ return exp(0.1*(x+y)); }
-double jacobian( double x, double y) 
+double jacobian( double x, double y)
 {
-    return exp( x-M_PI)*(sin(x)+cos(x))*sin(y) * exp(y-M_PI)*sin(x)*(sin(y) + cos(y)) - sin(x)*exp(x-M_PI)*cos(y) * cos(x)*sin(y)*exp(y-M_PI); 
+    return exp( x-M_PI)*(sin(x)+cos(x))*sin(y) * exp(y-M_PI)*sin(x)*(sin(y) + cos(y)) - sin(x)*exp(x-M_PI)*cos(y) * cos(x)*sin(y)*exp(y-M_PI);
 }
 */
 
 dg::bc bcx = dg::PER;
 dg::bc bcy = dg::PER;
 double left( double x, double y) {return sin(x)*cos(y);}
 double right( double x, double y) {return cos(x)*sin(y);}
-double jacobian( double x, double y) 
+double jacobian( double x, double y)
 {
-    return cos(x)*cos(y)*cos(x)*cos(y) - sin(x)*sin(y)*sin(x)*sin(y); 
+    return cos(x)*cos(y)*cos(x)*cos(y) - sin(x)*sin(y)*sin(x)*sin(y);
 }
 ////These are for comparing to FD arakawa results
 //double left( double x, double y) {return sin(2.*M_PI*(x-hx/2.));}
@@ -60,7 +60,7 @@ int main()
 
     dg::ArakawaX<dg::CartesianGrid2d, Matrix, Vector> arakawa( grid);
     unsigned multi=20;
-    t.tic(); 
+    t.tic();
     for( unsigned i=0; i<multi; i++)
         arakawa( lhs, rhs, jac);
     t.toc();