diff --git a/hysop/operator/base/advection_dir.py b/hysop/operator/base/advection_dir.py
index addedf7b3a77340c66848950f1638893f7487e4d..ec8ffeabeadf6d5be5fbe0e6dd136cf3e01017db 100644
--- a/hysop/operator/base/advection_dir.py
+++ b/hysop/operator/base/advection_dir.py
@@ -1,5 +1,5 @@
from hysop.tools.numpywrappers import npw
-from hysop.tools.decorators import debug
+from hysop.tools.decorators import debug
from hysop.tools.types import check_instance, InstanceOf
from hysop.methods import Interpolation, TimeIntegrator, Remesh
from hysop.fields.continuous_field import Field, ScalarField
@@ -12,13 +12,14 @@ from hysop.core.memory.memory_request import MemoryRequest
from hysop.core.graph.computational_node import ComputationalGraphNode
from hysop.parameters.scalar_parameter import ScalarParameter
+
class DirectionalAdvectionBase(object):
__default_method = {
- TimeIntegrator: Euler,
- Interpolation: Interpolation.LINEAR,
- Remesh: Remesh.L2_1,
- }
+ TimeIntegrator: Euler,
+ Interpolation: Interpolation.LINEAR,
+ Remesh: Remesh.L2_1,
+ }
__available_methods = {
TimeIntegrator: InstanceOf(ExplicitRungeKutta),
@@ -31,6 +32,7 @@ class DirectionalAdvectionBase(object):
dm = super(DirectionalAdvectionBase, cls).default_method()
dm.update(cls.__default_method)
return dm
+
@classmethod
def available_methods(cls):
am = super(DirectionalAdvectionBase, cls).available_methods()
@@ -39,21 +41,21 @@ class DirectionalAdvectionBase(object):
@debug
def __new__(cls, velocity, relative_velocity, splitting_direction,
- advected_fields_in, advected_fields_out,
- variables, dt, velocity_cfl,
- remesh_criteria_eps=None,
- **kwds):
+ advected_fields_in, advected_fields_out,
+ variables, dt, velocity_cfl,
+ remesh_criteria_eps=None,
+ **kwds):
return super(DirectionalAdvectionBase, cls).__new__(cls,
- input_fields=None, output_fields=None,
- input_params=None, output_params=None,
- splitting_direction=splitting_direction, **kwds)
+ input_fields=None, output_fields=None,
+ input_params=None, output_params=None,
+ splitting_direction=splitting_direction, **kwds)
@debug
def __init__(self, velocity, relative_velocity, splitting_direction,
- advected_fields_in, advected_fields_out,
- variables, dt, velocity_cfl,
- remesh_criteria_eps=None,
- **kwds):
+ advected_fields_in, advected_fields_out,
+ variables, dt, velocity_cfl,
+ remesh_criteria_eps=None,
+ **kwds):
"""
Particle advection of field(s) in a given direction,
on any backend, with cartesian remeshing.
@@ -75,8 +77,8 @@ class DirectionalAdvectionBase(object):
Dictionary of continuous fields as keys and topologies as values.
dt: ScalarParameter
Timestep parameter that will be used for time integration.
- velocity_cfl: float
- Velocity cfl in given direction.
+ velocity_cfl: float or tuple
+ Velocity cfl in given direction. In case of tuple given, value is taken from splitting direction.
remesh_criteria_eps: int
Minimum number of epsilons (in specified precision) that will trigger
remeshing. By default every non zero value is remeshed.
@@ -96,33 +98,36 @@ class DirectionalAdvectionBase(object):
check_instance(advected_fields_in, tuple, values=Field)
check_instance(advected_fields_out, tuple, values=Field)
check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
- check_instance(velocity_cfl, (float,int))
+ check_instance(velocity_cfl, (float, int, tuple))
+ if isinstance(velocity_cfl, tuple):
+ velocity_cfl = velocity_cfl[velocity.domain.dim-1-splitting_direction]
+ check_instance(velocity_cfl, (float, int))
check_instance(dt, ScalarParameter)
- check_instance(relative_velocity, tuple, values=(str,float), size=velocity.nb_components)
+ check_instance(relative_velocity, tuple, values=(str, float), size=velocity.nb_components)
velocity_cfl = float(velocity_cfl)
- assert (len(advected_fields_in)==len(advected_fields_out)), '|inputs| != |outputs|'
- assert (velocity_cfl>0.0), 'velocity_cfl cfl <= 0'
+ assert (len(advected_fields_in) == len(advected_fields_out)), '|inputs| != |outputs|'
+ assert (velocity_cfl > 0.0), 'velocity_cfl cfl <= 0'
# expand tensors
advected_scalars_in, advected_scalars_out = (), ()
input_tensors, output_tensors = (), ()
for (Sin, Sout) in zip(advected_fields_in, advected_fields_out):
- msg='Input and output field shape mismatch, got field {} of shape {} '
- msg+='and field {} of shape {}.'
+ msg = 'Input and output field shape mismatch, got field {} of shape {} '
+ msg += 'and field {} of shape {}.'
if (Sin.is_tensor ^ Sout.is_tensor):
if Sin.is_tensor:
- msg=msg.format(Sin.short_description(), Sin.shape,
- Sout.short_description(), '(1,)')
+ msg = msg.format(Sin.short_description(), Sin.shape,
+ Sout.short_description(), '(1,)')
else:
- msg=msg.format(Sin.short_description(), '(1,)',
- Sout.short_description(), Sout.shape)
+ msg = msg.format(Sin.short_description(), '(1,)',
+ Sout.short_description(), Sout.shape)
raise RuntimeError(msg)
if (Sin.is_tensor and Sout.is_tensor):
if (Sin.shape != Sout.shape):
- msg=msg.format(Sin.short_description(), Sin.shape,
- Sout.short_description(), Sout.shape)
+ msg = msg.format(Sin.short_description(), Sin.shape,
+ Sout.short_description(), Sout.shape)
raise RuntimeError(msg)
input_tensors += (Sin,)
output_tensors += (Sout,)
@@ -131,8 +136,8 @@ class DirectionalAdvectionBase(object):
check_instance(advected_scalars_in, tuple, values=ScalarField)
check_instance(advected_scalars_out, tuple, values=ScalarField)
- assert len(set(advected_scalars_in))==len(advected_scalars_in), \
- '|sinputs| != |soutputs|'
+ assert len(set(advected_scalars_in)) == len(advected_scalars_in), \
+ '|sinputs| != |soutputs|'
is_inplace = True
for (ifield, ofield) in zip(advected_scalars_in, advected_scalars_out):
@@ -144,31 +149,31 @@ class DirectionalAdvectionBase(object):
if velocity.is_tensor:
Vd = velocity[splitting_direction]
else:
- assert splitting_direction==0
+ assert splitting_direction == 0
Vd = velocity
Ud = relative_velocity[splitting_direction]
- input_fields = { Vd: variables[velocity] }
+ input_fields = {Vd: variables[velocity]}
output_fields = {}
- input_params = { dt.name: dt }
+ input_params = {dt.name: dt}
output_params = {}
for (ifield, ofield) in zip(advected_fields_in, advected_fields_out):
- input_fields[ifield] = variables[ifield]
+ input_fields[ifield] = variables[ifield]
output_fields[ofield] = variables[ofield]
super(DirectionalAdvectionBase, self).__init__(
- input_fields=input_fields,
- output_fields=output_fields,
- input_params=input_params,
- output_params=output_params,
- splitting_direction=splitting_direction,
- **kwds)
+ input_fields=input_fields,
+ output_fields=output_fields,
+ input_params=input_params,
+ output_params=output_params,
+ splitting_direction=splitting_direction,
+ **kwds)
self.velocity = Vd
self.relative_velocity = Ud
- self.first_scalar = advected_scalars_in[0]
- self.advected_fields_in = advected_scalars_in
+ self.first_scalar = advected_scalars_in[0]
+ self.advected_fields_in = advected_scalars_in
self.advected_fields_out = advected_scalars_out
self.dt = dt
@@ -176,33 +181,33 @@ class DirectionalAdvectionBase(object):
self.is_inplace = is_inplace
@debug
- def handle_method(self,method):
+ def handle_method(self, method):
super(DirectionalAdvectionBase, self).handle_method(method)
remesh_kernel = method.pop(Remesh)
if isinstance(remesh_kernel, Remesh):
remesh_kernel = RemeshKernel.from_enum(remesh_kernel)
- self.remesh_kernel = remesh_kernel
- self.interp = method.pop(Interpolation)
+ self.remesh_kernel = remesh_kernel
+ self.interp = method.pop(Interpolation)
self.time_integrator = method.pop(TimeIntegrator)
@debug
def get_field_requirements(self):
requirements = super(DirectionalAdvectionBase, self).get_field_requirements()
- direction = self.splitting_direction
+ direction = self.splitting_direction
is_inplace = self.is_inplace
- velocity = self.velocity
- velocity_cfl = self.velocity_cfl
+ velocity = self.velocity
+ velocity_cfl = self.velocity_cfl
v_topo, v_requirements = requirements.get_input_requirement(velocity)
try:
- v_dx = v_topo.space_step
+ v_dx = v_topo.space_step
except AttributeError:
- v_dx = v_topo.mesh.space_step
+ v_dx = v_topo.mesh.space_step
- dim = velocity.domain.dim
+ dim = velocity.domain.dim
ndir = dim-1-direction
scalar = self.first_scalar
@@ -267,8 +272,8 @@ class DirectionalAdvectionBase(object):
def scalars_out_cache_ghosts(cls, scalar_cfl, remesh_kernel):
"""Return the minimum number of ghosts for remeshed scalars."""
assert scalar_cfl > 0.0, 'cfl <= 0.0'
- assert remesh_kernel.n >=1 , 'Bad remeshing kernel.'
- if remesh_kernel.n >1:
+ assert remesh_kernel.n >= 1, 'Bad remeshing kernel.'
+ if remesh_kernel.n > 1:
assert remesh_kernel.n % 2 == 0, 'Odd remeshing kernel moments.'
min_ghosts = int(npw.floor(scalar_cfl)+1+remesh_kernel.n//2)
return min_ghosts
@@ -279,26 +284,26 @@ class DirectionalAdvectionBase(object):
dvelocity = self.input_discrete_fields[self.velocity]
# scalar field -> discrete_field
- dadvected_fields_in = {ifield: self.input_discrete_fields[ifield]
- for ifield in self.advected_fields_in}
+ dadvected_fields_in = {ifield: self.input_discrete_fields[ifield]
+ for ifield in self.advected_fields_in}
dadvected_fields_out = {ofield: self.output_discrete_fields[ofield]
- for ofield in self.advected_fields_out}
+ for ofield in self.advected_fields_out}
- self.dvelocity = dvelocity
- self.dadvected_fields_in = dadvected_fields_in
+ self.dvelocity = dvelocity
+ self.dadvected_fields_in = dadvected_fields_in
self.dadvected_fields_out = dadvected_fields_out
self.is_bilevel = None
- if any(self.dvelocity.compute_resolution != \
+ if any(self.dvelocity.compute_resolution !=
next(iter(self.dadvected_fields_out.values())).compute_resolution):
self.is_bilevel = self.dvelocity.compute_resolution
@debug
def get_work_properties(self):
- requests = super(DirectionalAdvectionBase, self).get_work_properties()
+ requests = super(DirectionalAdvectionBase, self).get_work_properties()
f = next(iter(self.dadvected_fields_in.values()))
(pos, request, request_id) = MemoryRequest.cartesian_dfield_like(name='position', dfield=f,
- ghosts=0, nb_components=1, is_read_only=False)
+ ghosts=0, nb_components=1, is_read_only=False)
requests.push_mem_request(request_id, request)
assert all(f.compute_resolution == pos.resolution)
self.dposition = pos
@@ -318,7 +323,7 @@ class DirectionalAdvectionBase(object):
raise ValueError('work is None.')
self.dposition.honor_memory_request(work, self)
- ## Backend methods
+ # Backend methods
# DirectionalOperatorBase
@classmethod
def supported_dimensions(cls):