GPU rolling sliding solver added

Author: mikeWShef

examples/Quasi Static Steps - Normal contact with movement.ipynb

@@ -1,5 +1,6 @@
 import multiprocessing
 import os
+import numpy
 """Top-level package for SlipPY."""
 
 __author__ = """Michael Watson"""
@@ -9,19 +10,25 @@
 try:
     import cupy  # noqa: F401
     CUDA = True
-    asnumpy = cupy.asnumpy
     xp = cupy
 except ImportError:
     CUDA = False
-    import numpy
-    asnumpy = numpy.asarray
     xp = numpy
 
+
+def asnumpy(obj):
+    if CUDA:
+        return cupy.asnumpy(obj)
+    return numpy.asarray(obj)
+
+
 CORES = multiprocessing.cpu_count()
 OUTPUT_DIR = os.getcwd()
 ERROR_IF_MISSING_MODEL = True
 ERROR_IF_MISSING_SUB_MODEL = True
 ERROR_IN_DATA_CHECK = True
+dtype = 'float64'
+material_names = []
 
 
 class OverRideCuda:

examples/Recreating the hertz solution numerically (normal contact).ipynb

@@ -116,7 +116,7 @@
 from .lubrication_steps import IterSemiSystem
 from .models import ContactModel
 from .static_step import StaticStep
-# from .steps import InitialStep
+from .steps import _ModelStep, RepeatingStateStep
 from .unified_reynolds_solver import UnifiedReynoldsSolver
 from .quasi_static_step import QuasiStaticStep
 from . import sub_models
@@ -125,5 +125,5 @@
            'lubricant_models', 'IterSemiSystem', 'Elastic', 'Rigid', 'rigid', 'elastic_influence_matrix',
            'ContactModel', 'OutputRequest', 'OutputReader', 'OutputSaver', 'read_output',
            'StaticStep', 'UnifiedReynoldsSolver', 'sub_models', 'QuasiStaticStep', 'sub_models',
-           'guess_loads_from_displacement', 'bccg', 'plan_convolve', 'plan_multi_convolve'
-           ]
+           'guess_loads_from_displacement', 'bccg', 'plan_convolve', 'plan_multi_convolve', '_ModelStep',
+           'RepeatingStateStep']

slippy/__init__.py

@@ -646,7 +646,6 @@ def hertz_full(r1: typing.Union[typing.Sequence, float], r2: typing.Union[typing
         zeta = [0 if v1 <= 0.1938 else 0.223 + 2.321 * v1 - 2.397 * v1 ** 2 for v1 in v]
         results['max_von_mises_stress'] = [c1 * p0 for c1 in c]
         results['max_von_mises_depth'] = [a * zeta1 for zeta1 in zeta]
-        results['total_deflection'] = 1 / np.pi / e_star * load * (np.log(4 * np.pi * e_star * r / load) - 1)
 
         # The following is taken from deeg
 

slippy/contact/__init__.py

@@ -85,6 +85,9 @@ class IterSemiSystem(_ModelStep):
         periodic in that direction. NOTE: this only controls the deformation result from the materials, ensure that the
         reynolds equation solver used supports periodic solutions and is set to produce a periodic solution to the
         pressure equation. Additionally, ensure that the axes of periodicity match.
+    periodic_im_repeats: tuple, optional (1,1)
+        The number of times the influence matrix should be wrapped along periodic dimensions, only used if at least one
+        of periodic axes is True. This is necessary to ensure truly periodic behaviour, no physical limit exists
     max_it_pressure: int, optional (100)
         The maximum number of iterations in the fluid pressure calculation loop
     rtol_pressure: float, optional (1e-7)
@@ -203,6 +206,7 @@ def __init__(self, step_name: str, reynolds_solver: _NonDimensionalReynoldSolver
                  movement_interpolation_mode: str = 'linear',
                  profile_interpolation_mode: str = 'nearest',
                  periodic_geometry: bool = False, periodic_axes: tuple = (False, False),
+                 periodic_im_repeats: tuple = (1, 1),
                  max_it_pressure: int = 5000, rtol_pressure: float = 2e-6,
                  max_it_interference: int = 5000,
                  rtol_interference: float = 1e-4,
@@ -218,7 +222,7 @@ def __init__(self, step_name: str, reynolds_solver: _NonDimensionalReynoldSolver
         self.profile_interpolation_mode = profile_interpolation_mode
         self._periodic_profile = periodic_geometry
         self._periodic_axes = periodic_axes
-
+        self._periodic_im_repeats = periodic_im_repeats
         self._max_it_pressure = max_it_pressure
         self._max_it_interference = max_it_interference
         self._rtol_pressure = rtol_pressure
@@ -367,7 +371,7 @@ def solve(self, previous_state: dict, output_file):
         for s in self.sub_models:
             s.no_time = self._no_time
 
-        relative_time = np.linspace(0, 1, self.number_of_steps)
+        relative_time = np.linspace(0, 1, self.number_of_steps+1)[1:]
         just_touching_gap = None
 
         original = dict()
@@ -396,14 +400,14 @@ def solve(self, previous_state: dict, output_file):
 
             # make a new loads function if we need it
             if (previous_gap_shape is None or previous_gap_shape != just_touching_gap.shape) and im_mats:
-                span = just_touching_gap.shape
+                span = tuple([s*(2-pa) for s, pa in zip(just_touching_gap.shape, self._periodic_axes)])
                 max_pressure = self.reynolds.dimensionalise_pressure(min([surf_1_material.max_load,
                                                                           surf_2_material.max_load]), True)
                 self._nd_max_pressure = max_pressure
-                im1 = surf_1_material.influence_matrix(span=span, grid_spacing=[gs] * 2,
-                                                       components=['zz'])['zz']
-                im2 = surf_2_material.influence_matrix(span=span, grid_spacing=[gs] * 2,
-                                                       components=['zz'])['zz']
+                im1 = surf_1_material.influence_matrix(components=['zz'], grid_spacing=[gs] * 2,
+                                                       span=span, periodic_strides=self._periodic_im_repeats)['zz']
+                im2 = surf_2_material.influence_matrix(components=['zz'], grid_spacing=[gs] * 2, span=span,
+                                                       periodic_strides=self._periodic_im_repeats)['zz']
                 total_im = im1 + im2
                 loads_func = plan_convolve(just_touching_gap, total_im, circular=self._periodic_axes)
                 previous_gap_shape = just_touching_gap.shape

slippy/contact/_step_utils.py

@@ -71,6 +71,8 @@ def __init__(self, name: str, surface_1: _SurfaceABC, surface_2: _SurfaceABC = N
         if lubricant is not None:
             self.lubricant_model = lubricant
         self.steps = OrderedDict({'Initial': InitialStep()})
+        self.current_step = None
+        self.current_step_start_time = None
         if output_dir is not None:
             if not os.path.isdir(output_dir):
                 try:
@@ -207,7 +209,7 @@ def solve(self, verbose: bool = False, skip_data_check: bool = False):
         if os.path.exists(self.log_file_name):
             os.remove(self.log_file_name)
 
-        current_state = None
+        current_state = {'time': 0.0}
 
         with ExitStack() as stack:
             output_writer = stack.enter_context(OutputSaver(self.name))
@@ -222,6 +224,8 @@ def solve(self, verbose: bool = False, skip_data_check: bool = False):
 
             for this_step in self.steps:
                 print(f"Solving step {this_step}")
+                self.current_step = self.steps[this_step]
+                self.current_step_start_time = current_state['time']
                 for output in self.steps[this_step].outputs:
                     output.new_step(current_state['time'])
                 current_state = self.steps[this_step].solve(current_state, output_writer)

slippy/contact/hertz.py

@@ -3,8 +3,6 @@
 import typing
 import warnings
 
-from scipy.optimize import root_scalar
-
 from .steps import _ModelStep
 from ._model_utils import get_gap_from_model
 from ._step_utils import HeightOptimisationFunction, make_interpolation_func
@@ -81,6 +79,14 @@ class QuasiStaticStep(_ModelStep):
     periodic_axes: tuple, optional ((False, False))
         For each True value the corresponding axis will be solved by circular convolution, meaning the result is
         periodic in that direction
+    periodic_im_repeats: tuple, optional (1,1)
+        The number of times the influence matrix should be wrapped along periodic dimensions, only used if at least one
+        of periodic axes is True. This is necessary to ensure truly periodic behaviour, no physical limit exists
+    method: {'auto', 'pk', 'double'}, optional ('auto')
+        The method by which the normal contact is solved, only used for load controlled contact.
+        'pk' uses the Polonsky and Keer algorithm for elastic contact.
+        'double' uses a double iteration procedure, suitable for elastic contact with a maximum pressure.
+        'auto' automatically selects 'pk' if there is no maximum pressure and 'double' if there is.
     max_it_interference: int, optional (100)
         The maximum number of iterations used to find the interference between the surfaces, only used if
         a total normal load is specified (Not used if contact is displacement controlled)
@@ -156,6 +162,7 @@ def __init__(self, step_name: str, number_of_steps: int, no_time: bool = False,
                  movement_interpolation_mode: str = 'linear',
                  profile_interpolation_mode: str = 'nearest',
                  periodic_geometry: bool = False, periodic_axes: tuple = (False, False),
+                 periodic_im_repeats: tuple = (1, 1), method: str = 'auto',
                  max_it_interference: int = 100, rtol_interference=1e-3,
                  max_it_displacement: int = None, rtol_displacement=1e-5, no_update_warning: bool = True,
                  upper: float = 4.0):
@@ -168,6 +175,7 @@ def __init__(self, step_name: str, number_of_steps: int, no_time: bool = False,
         self.total_time = time_period
         if not no_time and fast_ld:
             raise ValueError("Cannot have time dependence and fast_ld, either set no_time True or set fast_ld False")
+        self._periodic_im_repeats = periodic_im_repeats
         self._fast_ld = fast_ld
         self._relative_loading = relative_loading
         self.profile_interpolation_mode = profile_interpolation_mode
@@ -178,6 +186,11 @@ def __init__(self, step_name: str, number_of_steps: int, no_time: bool = False,
         self._max_it_displacement = max_it_displacement
         self._rtol_displacement = rtol_displacement
         self.number_of_steps = number_of_steps
+
+        if method not in {'auto', 'pk', 'double'}:
+            raise ValueError(f"Unrecognised method for step {step_name}: {method}")
+
+        self._method = method
         self._height_optimisation_func = None
         self._adhesion = adhesion
         self._unloading = unloading
@@ -256,7 +269,7 @@ def solve(self, previous_state, output_file):
         else:  # displacement controlled
             update_func = self._solve_displacement_controlled
 
-        relative_time = np.linspace(0, 1, self.number_of_steps)
+        relative_time = np.linspace(0, 1, self.number_of_steps + 1)[1:]
         just_touching_gap = None
 
         original = dict()
@@ -314,7 +327,7 @@ def solve(self, previous_state, output_file):
     @property
     def upper(self):
         if self._upper is None:
-            self._upper = np.max(self._just_touching_gap)*self._upper_factor
+            self._upper = np.max(self._just_touching_gap) * self._upper_factor
         return self._upper
 
     def update_movement(self, relative_time, original):
@@ -327,47 +340,53 @@ def update_movement(self, relative_time, original):
     def _solve_load_controlled(self, current_state) -> dict:
         # if there is time dependence or we don't already have one, make a new height optimiser
         if not self._no_time or self._height_optimisation_func is None:
-            h_opt_func = HeightOptimisationFunction(just_touching_gap=self._just_touching_gap,
-                                                    model=self.model,
-                                                    adhesion_model=self._adhesion_model,
-                                                    initial_contact_nodes=self._initial_contact_nodes,
-                                                    max_it_inner=self._max_it_displacement,
-                                                    tol_inner=self._rtol_displacement,
-                                                    material_options=dict(),
-                                                    max_set_load=self.normal_load,
-                                                    tolerance=self._rtol_interference,
-                                                    periodic_axes=self._periodic_axes)
-            self._height_optimisation_func = h_opt_func
+            opt_func = HeightOptimisationFunction(just_touching_gap=self._just_touching_gap,
+                                                  model=self.model,
+                                                  adhesion_model=self._adhesion_model,
+                                                  initial_contact_nodes=self._initial_contact_nodes,
+                                                  max_it_inner=self._max_it_displacement,
+                                                  tol_inner=self._rtol_displacement,
+                                                  material_options=dict(),
+                                                  max_set_load=self.normal_load,
+                                                  tolerance=self._rtol_interference,
+                                                  periodic_axes=self._periodic_axes,
+                                                  periodic_im_repeats=self._periodic_im_repeats)
+            self._height_optimisation_func = opt_func
         else:
-            h_opt_func = self._height_optimisation_func
+            opt_func = self._height_optimisation_func
 
         if self._unloading and 'contact_nodes' in current_state:
             contact_nodes = current_state['contact_nodes']
         else:
             contact_nodes = None
             # contact_nodes = np.ones(self._just_touching_gap.shape, dtype=np.bool)
+        if self._method == 'auto':
+            if np.isinf(opt_func.max_pressure):
+                self._method = 'pk'
+            else:
+                self._method = 'double'
 
-        h_opt_func.change_load(self.normal_load, contact_nodes)
-
-        # need to set bounds and pick a sensible starting point
-        upper = self.upper
-        print(f'upper bound set at: {upper}')
-        if self._no_time:
-            brackets = h_opt_func.get_bounds_from_cache(0, upper)
+        if self._method == 'pk':
+            opt_func.contact_nodes = None
+            opt_func.p_and_k(self.normal_load)
         else:
-            brackets = (0, upper)
-        print(f'Bounds adjusted using cache to: {brackets}')
-        print(f'Interference tolerance set to {self._rtol_interference} Relative')
+            opt_func.change_load(self.normal_load, contact_nodes)
+            # need to set bounds and pick a sensible starting point
+            upper = self.upper
+            print(f'upper bound set at: {upper}')
+            if self._no_time:
+                brackets = opt_func.get_bounds_from_cache(0, upper)
+            else:
+                brackets = (0, upper)
+            print(f'Bounds adjusted using cache to: {brackets}')
+            print(f'Interference tolerance set to {self._rtol_interference} Relative')
+            opt_func.brent(0, upper, r_tol=self._rtol_interference, max_iter=self._max_it_interference)
 
-        opt_result = root_scalar(h_opt_func, bracket=brackets, rtol=self._rtol_interference,
-                                 maxiter=self._max_it_interference, args=(current_state,))
         # noinspection PyProtectedMember
-        if h_opt_func._results is None:
-            h_opt_func((brackets[0]+brackets[1])/2, current_state)
-        results = h_opt_func.results
-        results['interference'] = opt_result.root
-        load_conv = (np.abs(results['total_normal_load']-self.normal_load) / self.normal_load) < 0.05
-        results['converged'] = bool(load_conv) and not h_opt_func.last_call_failed
+
+        results = opt_func.results
+        load_conv = (np.abs(results['total_normal_load'] - self.normal_load) / self.normal_load) < 0.05
+        results['converged'] = bool(load_conv) and not opt_func.last_call_failed
         return results
 
     def _solve_displacement_controlled(self, current_state):
@@ -381,7 +400,8 @@ def _solve_displacement_controlled(self, current_state):
                                                     material_options=dict(),
                                                     max_set_load=1,
                                                     tolerance=self._rtol_interference,
-                                                    periodic_axes=self._periodic_axes)
+                                                    periodic_axes=self._periodic_axes,
+                                                    periodic_im_repeats=self._periodic_im_repeats)
             self._height_optimisation_func = h_opt_func
         else:
             h_opt_func = self._height_optimisation_func