Add SCALES (lib) and multiscale tests for advection

ci/scripts/config.sh

@@ -9,12 +9,12 @@ fi
 if [ -d "$1" ]; then
     echo "Folder $1 already exists."
     exit 1
-fi 
+fi
 
 if [ -d "$2" ]; then
     echo "Folder $2 already exists."
     exit 1
-fi 
+fi
 
 ROOT_DIR="$(pwd)"
 BUILD_DIR="$1"
@@ -22,7 +22,7 @@ INSTALL_DIR="$2"
 
 mkdir -p $BUILD_DIR
 cd $BUILD_DIR
-CC="$3" CXX="$4" FC="$5" cmake -DCMAKE_BUILD_TYPE=Release -DVERBOSE=OFF -DHYSOP_INSTALL=$INSTALL_DIR $ROOT_DIR
+CC="$3" CXX="$4" FC="$5" cmake -DCMAKE_BUILD_TYPE=Release -DVERBOSE=OFF -DWITH_SCALES=ON -DHYSOP_INSTALL=$INSTALL_DIR $ROOT_DIR
 
 if [ ! -f Makefile ]; then
     echo "The makefile has not been generated."

ci/scripts/test.sh

@@ -65,6 +65,10 @@ python "$HYSOP_DIR/operator/tests/test_directional_diffusion.py"
 python "$HYSOP_DIR/operator/tests/test_diffusion.py"
 python "$HYSOP_DIR/operator/tests/test_directional_stretching.py"
 
+# If scales (fortran advection library) is installed
+python -c "from hysop.f2hysop import scales2py as scales" && python "$HYSOP_DIR/operator/tests/test_multiscale_advection.py"
+python -c "from hysop.f2hysop import scales2py as scales" && python "$HYSOP_DIR/operator/tests/test_scales_advection.py"
+
 if [ "$HAS_CACHE_DIR" = true ]; then
     cp -r /root/.cache/* $CACHE_DIR/
     find $CACHE_DIR -name '*.lock' -delete

hysop/backend/host/fortran/operator/scales_advection.py

@@ -255,12 +255,16 @@ class ScalesAdvection(ComputationalGraphOperator):
                     # (slower than basic): solve_advection_inter
                 else:
                     self._scales_func.append(scales.solve_advection_vect)
-            else:
+            elif df.nb_components == 1:
                 if self.is_multiscale is not None:
                     self._scales_func.append(
                         scales.solve_advection_inter_basic)
                 else:
                     self._scales_func.append(scales.solve_advection)
+            else:
+                msg = "SCALES library for advection works only with "
+                msg += "3D scalar fields or 3D 3-components vectors."
+                raise RuntimeError(msg)
 
     ## Backend methods
     # DirectionalOperatorBase

hysop/operator/tests/test_scales_advection.py

@@ -5,13 +5,13 @@ from hysop.testsenv import iter_clenv
 from hysop.tools.contexts import printoptions
 from hysop.tools.types import first_not_None, to_tuple
 from hysop.tools.io_utils import IO
-from hysop.operator.advection  import Advection  # Scales fortran advection
 from hysop.operator.directional.advection_dir import DirectionalAdvection
+from hysop.operator.advection import Advection
 from hysop.parameters.scalar_parameter import ScalarParameter
 
-from hysop import Field, Box, Problem
+from hysop import Field, Box
 from hysop.methods import Remesh, TimeIntegrator
-from hysop.constants import Implementation, DirectionLabels
+from hysop.constants import Implementation, DirectionLabels, HYSOP_REAL
 from hysop.numerics.splitting.strang import StrangSplitting, StrangOrder
 from hysop.numerics.odesolvers.runge_kutta import Euler, RK2, RK4
 from hysop.numerics.remesh.remesh import RemeshKernel
@@ -51,19 +51,25 @@ class TestDirectionalAdvectionOperator(object):
 
         shape = tuple(npw.random.randint(low=size_min, high=size_max, size=dim).tolist())
 
-        flt_types = (npw.float64,)
-        time_integrators = (RK2,)
-        remesh_kernels = (Remesh.L2_1,)
-        velocity_cfls = (1.89,)
+        if self.enable_extra_tests:
+            flt_types = (HYSOP_REAL, )
+            time_integrators = (RK2, )
+            remesh_kernels =  (Remesh.L2_1, Remesh.L4_2, Remesh.L6_4)
+            velocity_cfls = (0.62, 1.89)
+        else:
+            flt_types = (HYSOP_REAL,)
+            time_integrators = (RK2,)
+            remesh_kernels = (Remesh.L4_2,)
+            velocity_cfls = (1.89,)
 
         domain = Box(length=(2*npw.pi,)*dim)
         for dtype in flt_types:
             Vin  = Field(domain=domain, name='Vin', dtype=dtype,
                     nb_components=dim, register_object=False)
             Sin  = Field(domain=domain, name='Sin', dtype=dtype,
-                    nb_components=1, register_object=False)
+                    nb_components=3, register_object=False)
             Sout = Field(domain=domain, name='Sout', dtype=dtype,
-                    nb_components=1, register_object=False)
+                    nb_components=3, register_object=False)
             for time_integrator in time_integrators:
                 for remesh_kernel in remesh_kernels:
                     for velocity_cfl in velocity_cfls:
@@ -112,16 +118,17 @@ class TestDirectionalAdvectionOperator(object):
         dt.value = (0.99 * velocity_cfl) / (max(shape)-1)
 
         ref_impl = Implementation.PYTHON
-        implementations = [ref_impl, Implementation.FORTRAN]
+        implementations = Advection.implementations().keys()
+        implementations = [ref_impl] + implementations
 
         method = {TimeIntegrator: time_integrator, Remesh: remesh_kernel}
 
         def iter_impl(impl):
-            print impl
+            base_kwds = dict(velocity=vin, advected_fields=sin, advected_fields_out=sout, dt=dt,
+                             variables=variables, implementation=impl,
+                             method=method, name='advection_{}'.format(str(impl).lower()))
             if impl is Implementation.PYTHON:
-                da = DirectionalAdvection(velocity=vin, advected_fields=sin, advected_fields_out=sout, dt=dt,
-                                          velocity_cfl=velocity_cfl, variables=variables, implementation=impl,
-                                          method=method, name='advection_{}'.format(str(impl).lower()))
+                da = DirectionalAdvection(velocity_cfl=velocity_cfl, **base_kwds)
                 split = StrangSplitting(splitting_dim=dim,
                                         order=StrangOrder.STRANG_SECOND_ORDER)
                 split.push_operators(da)
@@ -129,10 +136,9 @@ class TestDirectionalAdvectionOperator(object):
                     split.push_copy(dst=sin, src=sout)
                 yield 'default', split
             elif impl is Implementation.FORTRAN:
-                da = Advection(velocity=vin, advected_fields=sin, advected_fields_out=sout, dt=dt,
-                               variables=variables, implementation=impl,
-                               method=method, name='advection_{}'.format(str(impl).lower()))
-                yield 'Scales', da.to_graph()
+                assert dim==3, "Scales is only 3D"
+                da = Advection(**base_kwds)
+                yield 'SCALES', da.to_graph()
             else:
                 msg='Unknown implementation to test {}.'.format(impl)
                 raise NotImplementedError(msg)
@@ -176,7 +182,7 @@ class TestDirectionalAdvectionOperator(object):
                             if (k>0):
                                 op.apply()
 
-                            output = tuple(dsout[i][dsout.compute_slices].get().handle.copy()
+                            output = tuple(dsout[i].get().handle.copy()[dsout.compute_slices]
                                         for i in xrange(dsout.nb_components))
 
                             for i in xrange(dsout.nb_components):
@@ -208,7 +214,7 @@ class TestDirectionalAdvectionOperator(object):
                                 for i in xrange(dsout.nb_components):
                                     di = npw.abs(reference[i] - output[i])
                                     max_di = npw.max(di)
-                                    neps = 500
+                                    neps = 10000
                                     max_tol = neps*npw.finfo(dsout.dtype).eps
                                     if (max_di>max_tol):
                                         print '\nComparisson failed at step {} and component {}:'.format(k,i)
@@ -244,10 +250,9 @@ class TestDirectionalAdvectionOperator(object):
                         sys.stdout.flush()
                         raise
 
-
     def perform_tests(self):
+        #Scales is only 3D
         self._test(dim=3, is_inplace=True)
-
         print