Skip to content
Snippets Groups Projects
Commit 6bcec114 authored by EXT Jean-Matthieu Etancelin's avatar EXT Jean-Matthieu Etancelin
Browse files

add directional maximum cfl for advection

parent 569e6cab
No related branches found
No related tags found
2 merge requests!24Resolve "Add python3.x support",!15WIP: Resolve "HySoP with tasks"
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):
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment