added diffusion to poisson curl

examples/particles_above_salt/particles_above_salt_bc.py

@@ -48,7 +48,8 @@ def compute(args):
                                EnstrophyParameter, TimeParameters, \
                                VolumicIntegrationParameter
     from hysop.constants import Implementation, AdvectionCriteria, \
-                                BoxBoundaryCondition, BoundaryCondition
+                                BoxBoundaryCondition, BoundaryCondition, \
+                                Backend
 
     from hysop.operators import DirectionalAdvection, DirectionalStretching,       \
                                 Diffusion, ComputeMeanField,                       \
@@ -181,13 +182,6 @@ def compute(args):
         raise RuntimeError(msg)
     
     #> Diffusion of vorticity, S and C
-    diffuse_W = Diffusion(implementation=impl,
-             name='diffuse_{}'.format(vorti.name),
-             pretty_name=u'diff{}'.format(vorti.pretty_name.decode('utf-8')),
-             nu = nu_W,
-             Fin = vorti,
-             variables = {vorti: npts},
-             dt=dt, **extra_op_kwds)
     diffuse_S = Diffusion(implementation=impl,
              name='diffuse_S',
              pretty_name='diffS',
@@ -226,6 +220,7 @@ def compute(args):
     #> Poisson operator to recover the velocity from the vorticity
     poisson = PoissonCurl(name='poisson', velocity=velo, vorticity=vorti, 
                             variables={velo:npts, vorti: npts},
+                            diffusion=nu_W, dt=dt,
                             implementation=impl, **extra_op_kwds)
     
     #> Operator to compute the infinite norm of the velocity
@@ -241,10 +236,12 @@ def compute(args):
     io_params = IOParams(filename='fields', frequency=args.dump_freq)
     dump_fields = HDF_Writer(name='dump',
                              io_params=io_params,
+                             force_backend=Backend.OPENCL,
                              variables={velo: npts, 
                                         vorti: npts,
                                         C: npts, 
-                                        S: npts})
+                                        S: npts},
+                             **extra_op_kwds)
 
     #> Operator to compute and save mean fields
     axes = list(range(0, dim-1))
@@ -287,7 +284,7 @@ def compute(args):
 
     problem = Problem(method=method)
     problem.insert(poisson, 
-                   diffuse_W, diffuse_S, diffuse_C,
+                   diffuse_S, diffuse_C,
                    splitting, 
                    dump_fields,
                    compute_mean_fields,

examples/particles_above_salt/particles_above_salt_bc_3d.py

@@ -75,8 +75,8 @@ def compute(args):
     if (dim==2):
         (Sc, tau, Vp, Rs, Xo, Xn, N) = (0.70,  25, 0.04, 2.0, (-600,0), (600,750), (1537, 512))
     elif (dim==3):
-        (Sc, tau, Vp, Rs, Xo, Xn, N) = (30.00, 25, 0.04, 2.0, (-110,0,0), (65,100,100), (3073, 1024, 1024))
-        N = (257,64,64)
+        (Sc, tau, Vp, Rs, Xo, Xn, N) = (7.00, 25, 0.04, 2.0, (-110,0,0), (65,100,100), (3073, 1024, 1024))
+        N = (513,128,128)
     else:
         raise NotImplementedError
 
@@ -183,13 +183,6 @@ def compute(args):
         raise RuntimeError(msg)
     
     #> Diffusion of vorticity, S and C
-    diffuse_W = Diffusion(implementation=impl,
-             name='diffuse_{}'.format(vorti.name),
-             pretty_name=u'diff{}'.format(vorti.pretty_name.decode('utf-8')),
-             nu = nu_W,
-             Fin = vorti,
-             variables = {vorti: npts},
-             dt=dt, **extra_op_kwds)
     diffuse_S = Diffusion(implementation=impl,
              name='diffuse_S',
              pretty_name='diffS',
@@ -228,6 +221,7 @@ def compute(args):
     #> Poisson operator to recover the velocity from the vorticity
     poisson = PoissonCurl(name='poisson', velocity=velo, vorticity=vorti, 
                             variables={velo:npts, vorti: npts},
+                            diffusion=nu_W, dt=dt,
                             implementation=impl, **extra_op_kwds)
     
     #> Operator to compute the infinite norm of the velocity
@@ -244,7 +238,8 @@ def compute(args):
     dump_fields = HDF_Writer(name='dump',
                              io_params=io_params,
                              force_backend=Backend.OPENCL,
-                             variables={velo[0]: npts, 
+                             variables={vorti: npts,
+                                        velo: npts, 
                                         C: npts, 
                                         S: npts},
                              **extra_op_kwds)
@@ -280,7 +275,7 @@ def compute(args):
 
     problem = Problem(method=method)
     problem.insert(poisson, dump_fields,
-                   diffuse_W, diffuse_S, diffuse_C,
+                   diffuse_S, diffuse_C,
                    splitting, 
                    min_max_U, min_max_W, adapt_dt)
     problem.build()

examples/taylor_green/taylor_green.py

@@ -35,7 +35,8 @@ def init_vorticity(data, coords, component=None):
 def compute(args):
     from hysop import Box, Simulation, Problem, MPIParams
     from hysop.defaults import VelocityField, VorticityField, \
-                               EnstrophyParameter, TimeParameters
+                               EnstrophyParameter, TimeParameters, \
+                               ViscosityParameter
     from hysop.constants import Implementation, AdvectionCriteria, StretchingCriteria
 
     from hysop.operators import DirectionalAdvection, DirectionalStretchingDiffusion, \
@@ -86,6 +87,7 @@ def compute(args):
     velo  = VelocityField(domain=box, dtype=args.dtype)
     vorti = VorticityField(velocity=velo)
     enstrophy = EnstrophyParameter(dtype=args.dtype)
+    viscosity = ViscosityParameter(dtype=args.dtype, initial_value=(1.0/args.Re), const=True)
     
     ### Build the directional operators
     if (impl is Implementation.FORTRAN):
@@ -124,34 +126,17 @@ def compute(args):
                  vorticity = vorti,
                  variables = {velo: npts, vorti: npts},
                  dt=dt, **extra_op_kwds)
-    #> Directional diffusion
-    if (impl is Implementation.OPENCL):
-        diffuse = None
-        diffuse_dir = DirectionalDiffusion(implementation=impl,
-                 name='diffuse',
-                 fields = vorti,
-                 coeffs = (1.0/args.Re),
-                 variables = {vorti: npts},
-                 dt=dt, **extra_op_kwds)
-    else:
-        diffuse = Diffusion(
-                implementation=impl,
-                name='diffuse',
-                Fin = vorti,
-                viscosity = (1.0/args.Re),
-                variables = {vorti: npts},
-                dt=dt, **extra_op_kwds)
-        diffuse_dir = None
     
     ### Build standard operators
     #> Poisson operator to recover the velocity from the vorticity
     poisson = PoissonCurl(name='poisson', velocity=velo, vorticity=vorti, 
                             variables={velo:npts, vorti: npts}, 
                             projection=args.reprojection_frequency,
+                            diffusion=viscosity, dt=dt,
                             implementation=impl, **extra_op_kwds)
     #> We ask to dump the outputs of this operator
-    #poisson.dump_outputs(fields=(vorti,), frequency=args.dump_freq)
-    #poisson.dump_outputs(fields=(velo,),  frequency=args.dump_freq)
+    poisson.dump_outputs(fields=(vorti,), frequency=args.dump_freq)
+    poisson.dump_outputs(fields=(velo,),  frequency=args.dump_freq)
     
     #> Operator to compute the infinite norm of the velocity
     if (impl is Implementation.FORTRAN):
@@ -171,7 +156,7 @@ def compute(args):
     
     #> Directional splitting operator subgraph
     splitting = StrangSplitting(splitting_dim=dim, order=args.strang_order)
-    splitting.push_operators(advec_dir, stretch_dir, diffuse_dir, min_max_gradU)
+    splitting.push_operators(advec_dir, stretch_dir, min_max_gradU)
 
     ### Adaptive timestep operator
     adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True)
@@ -269,8 +254,7 @@ def compute(args):
             }
     )
     problem = Problem(method=method)
-    problem.insert(poisson, 
-                   advec, splitting, diffuse,
+    problem.insert(poisson, advec, splitting,
                    min_max_U, min_max_W, enstrophy_op, 
                    adapt_dt, plot)
     problem.build()
@@ -379,7 +363,7 @@ if __name__=='__main__':
                         tstart=0.0, tend=20.01, 
                         dt=1e-5, 
                         cfl=0.5, lcfl=0.125,
-                        dump_freq=100, dump_times=(),
+                        dump_freq=0, dump_times=(),
                         Re=1600.0)
 
     parser.run(compute)

hysop/backend/device/opencl/operator/poisson_curl.py

@@ -38,8 +38,8 @@ class OpenClPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClSymbolic):
         Win  = tuple(Ft.output_symbolic_array('{}in_hat'.format(Ft.field.var_name)) for Ft in W_Ft)
         Wout = tuple(Bt.input_symbolic_array('{}out_hat'.format(Bt.field.var_name)) if (Bt is not None) else None for Bt in W_Bt)
         Uout = tuple(Bt.input_symbolic_array('{}out_hat'.format(Bt.field.var_name)) for Bt in U_Bt)
-        K  = tuple(tuple(self.tg._indexed_wave_numbers[kd] for kd in kd1[::-1]) for kd1 in kd1s)
-        KK = tuple(tuple(self.tg._indexed_wave_numbers[kd] for kd in kd2[::-1]) for kd2 in kd2s)
+        K    = tuple(tuple(self.tg._indexed_wave_numbers[kd] for kd in kd1[::-1]) for kd1 in kd1s)
+        KK   = tuple(tuple(self.tg._indexed_wave_numbers[kd] for kd in kd2[::-1]) for kd2 in kd2s)
 
         def mul(Ki, other):
             if Ki.Wn.is_complex:
@@ -47,14 +47,21 @@ class OpenClPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClSymbolic):
             else:
                 return Ki*other
 
+        if self.should_diffuse:
+            nu, dt = self.nu.s, self.dt.s
+            for i,(win,wout,KKi) in enumerate(zip(Win,Wout,KK)):
+                F = sum(KKi)
+                expr = Assignment(wout, win / (1 - nu*dt*F))
+                self.require_symbolic_kernel('diffusion_kernel__{}'.format(i), expr)
+            Win = Wout
+
         exprs = ()
         indices = local_indices_symbols[:dim]
         cond = LogicalAND(*tuple(LogicalEQ(idx,0) for idx in indices))
         for i in xrange(wcomp):
-            e = sum(KK[i])
-            e = Assignment(Win[i], Select(Win[i]/e,0,cond))
-            exprs += (e,)
-        self.require_symbolic_kernel('poisson_kernel', *exprs)
+            F = sum(KK[i])
+            expr = Assignment(Win[i], Select(Win[i]/F,0,cond))
+            self.require_symbolic_kernel('poisson_kernel__{}'.format(i), expr)
         
         if (dim == 2):
             assert wcomp==1
@@ -81,34 +88,42 @@ class OpenClPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClSymbolic):
         self._build_projection_kernel()
         self._build_poisson_curl_kernel()
         self._build_ghost_exchangers()
+    
+    def _build_diffusion_kernel(self):
+        if self.should_diffuse:
+            diffusion_filters = ()
+            for i in xrange(self.W.nb_components):
+                knl, knl_kwds = \
+                        self.symbolic_kernels['diffusion_kernel__{}'.format(i)]
+                knl = functools.partial(knl, queue=self.cl_env.default_queue)
+                def F(knl=knl, knl_kwds=knl_kwds):
+                    return knl(**knl_kwds())
+                diffusion_filters += (F,)
+            self.diffusion_filters = diffusion_filters
+        else:
+            self.diffusion_filters = None
 
     def _build_poisson_curl_kernel(self):
-        filter_poisson, _ = self.symbolic_kernels['poisson_kernel']
-        self.poisson_filter = functools.partial(filter_poisson, queue=self.cl_env.default_queue)
-
-        F0, _ = self.symbolic_kernels['curl_kernel__0'] 
-        F1, _ = self.symbolic_kernels['curl_kernel__1'] 
-        if (self.dim==2):
-            curl_kernels = (F0,F1)
-        elif (self.dim==3):
-            F2, _ = self.symbolic_kernels['curl_kernel__2'] 
-            curl_kernels = (F0,F1,F2)
-        else:
-            msg='dim={}'.format(self.dim)
-            raise NotImplementedError(msg)
-        assert len(curl_kernels)==len(self.U_backward_transforms)
-        self.curl_filters = tuple(functools.partial(Fi, queue=self.cl_env.default_queue) for Fi in curl_kernels)
+        poisson_filters = ()
+        for i in xrange(self.W.nb_components):
+            knl, __ = self.symbolic_kernels['poisson_kernel__{}'.format(i)] 
+            Fi = functools.partial(knl, queue=self.cl_env.default_queue)
+            poisson_filters += (Fi,)
+        
+        curl_filters = ()
+        for i in xrange(self.U.nb_components):
+            knl, __ = self.symbolic_kernels['curl_kernel__{}'.format(i)] 
+            Fi = functools.partial(knl, queue=self.cl_env.default_queue)
+            curl_filters += (Fi,)
+        
+        self.poisson_filters = poisson_filters
+        self.curl_filters = curl_filters
 
     def _build_projection_kernel(self):
         if self.should_project:
+            raise NotImplementedError
             self.filter_projection, _ = self.symbolic_kernels['projection_kernel']
 
-    def _build_diffusion_kernel(self):
-        if self.should_diffuse:
-            filter_diffusion, diffusion_kernel_kwds = \
-                    self.symbolic_kernels['diffusion_kernel']
-            self.filter_diffusion = filter_diffusion
-
     def _build_ghost_exchangers(self):
         def noop():
             pass
@@ -132,21 +147,24 @@ class OpenClPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClSymbolic):
     def apply(self, simulation, **kwds):
         '''Solve the PoissonCurl equation.'''
         
-        diffuse=self.should_diffuse
-        project=self.do_project(simu=simulation)
+        diffuse = self.should_diffuse
+        project = self.do_project(simu=simulation)
         
         for Ft in self.W_forward_transforms:
             evt = Ft()
-        #if diffuse:
-            #evt = self.filter_diffusion(nu_dt)
-        #if project:
-            #evt = self.filter_projection()
-        #if (diffuse or project):
-            #for Bt in self.W_backward_transforms:
-                #evt = Bt()
-            #evt = self.exchange_W_ghosts()
-        evt = self.poisson_filter()
-        for (curl_filter,Bt) in zip(self.curl_filters, self.U_backward_transforms):
-            evt = curl_filter()
+        if diffuse:
+            for Fd in self.diffusion_filters:
+                evt = Fd()
+        if project:
+            evt = self.filter_projection()
+        if (diffuse or project):
+            for Bt in self.W_backward_transforms:
+                evt = Bt()
+            evt = self.exchange_W_ghosts()
+        for Fp in self.poisson_filters:
+            evt = Fp()
+        for (Fc, Bt) in zip(self.curl_filters, self.U_backward_transforms):
+            evt = Fc()
             evt = Bt()
         evt = self.exchange_U_ghosts()
+

hysop/operator/base/poisson_curl.py

@@ -60,9 +60,10 @@ class PoissonCurlOperatorBase(object):
 
         # check for diffusion
         should_diffuse = (diffusion is not None)
-        if should_diffuse and (dt is None):
-            msg='Diffusion has been specified but no timestep was given.'
-            raise RuntimeError(msg)
+        if should_diffuse:
+            if (dt is None):
+                msg='Diffusion has been specified but no timestep was given.'
+                raise RuntimeError(msg)
         else:
             dt = None
         
@@ -105,6 +106,7 @@ class PoissonCurlOperatorBase(object):
         input_fields  = { vorticity: wtopology }
         output_fields = { velocity:  vtopology }
         if should_diffuse:
+            assert (dt is not None), 'Diffusion timestep has not been given.'
             input_params[diffusion.name] = diffusion
             input_params[dt.name] = dt
         if (should_diffuse or should_project):
@@ -118,7 +120,7 @@ class PoissonCurlOperatorBase(object):
         self.dim = dim
         
         self.should_diffuse = should_diffuse
-        self.diffusion = diffusion
+        self.nu = diffusion
         self.dt = dt
 
         self.should_project = should_project