diff --git a/tests/mpi_routines.py b/tests/mpi_routines.py index 5a7fbda..91c73d2 100644 --- a/tests/mpi_routines.py +++ b/tests/mpi_routines.py @@ -1,308 +1,310 @@ # -*- coding: utf-8 -*- # # Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . -from __future__ import print_function import os import tempfile +import subprocess import tamaas as tm import numpy as np from mpi4py import MPI from numpy.linalg import norm from conftest import HertzFixture def make_surface(N): spectrum = tm.Isopowerlaw2D() spectrum.q0 = 2 spectrum.q1 = 2 spectrum.q2 = 16 spectrum.hurst = 0.8 generator = tm.SurfaceGeneratorRandomPhase2D(N) generator.spectrum = spectrum generator.random_seed = 1 return generator.buildSurface() def mpi_surface_generator(): N = [128, 128] tm.set_log_level(tm.LogLevel.debug) seq_surface = None comm = MPI.COMM_WORLD print('[{}] {}'.format(comm.rank, comm.size)) with tm.mpi.sequential(): if comm.rank == 0: seq_surface = make_surface(N) surface = make_surface(N) print('[{}] {}'.format(comm.rank, surface.shape)) recv = comm.gather(surface, root=0) if comm.rank == 0: gsurface = np.concatenate(recv) if False: import matplotlib.pyplot as plt plt.imshow(seq_surface) plt.colorbar() plt.figure() plt.imshow(gsurface) plt.colorbar() plt.show() # I think the only reason this assert works is because the # spectrum is compact in the Fourier domain -> surface is regular assert np.all(seq_surface == gsurface) def mpi_model_creation(): N = [20, 50, 50] S = [1., 1., 1.] comm = MPI.COMM_WORLD def get_discretizations(model): return model.shape, model.boundary_shape model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S[1:], N[1:]) n, bn = get_discretizations(model) n[0] = comm.allreduce(n[0]) bn[0] = comm.allreduce(bn[0]) assert n == N[1:] and bn == N[1:] model = tm.ModelFactory.createModel(tm.model_type.volume_2d, S, N) n, bn = get_discretizations(model) n[1] = comm.allreduce(n[1]) bn[0] = comm.allreduce(bn[0]) assert n == N and bn == N[1:] def mpi_polonsky_keer(): N = [1024, 1024] S = [1., 1.] load = 0.00001 comm = MPI.COMM_WORLD hertz = HertzFixture(N[0], load) surface = hertz.surface tm.set_log_level(tm.LogLevel.debug) def solve(): model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S, N) model.E, model.nu = 1, 0 local_surface = tm.mpi.scatter(surface) solver = tm.PolonskyKeerRey(model, local_surface, 1e-15) print('Created solver') solver.solve(load) return np.copy(model['traction']), np.copy(model['displacement']) tractions, displacements = map(tm.mpi.gather, solve()) if comm.rank == 0: perror = norm(tractions - hertz.pressure) / norm(hertz.pressure) derror = norm(displacements - hertz.displacement)\ / norm(hertz.displacement) if False: print(perror, derror) import matplotlib.pyplot as plt plt.imshow(tractions - hertz.pressure) plt.colorbar() plt.show() assert perror < 5e-3 assert derror < 8e-3 def mpi_polonsky_keer_compare(): N = [273, 273] S = [1., 1.] E, nu = 0.2, 0.3 load = 0.1 comm = MPI.COMM_WORLD seq_tractions = None rms = 0 with tm.mpi.sequential(): if comm.rank == 0: model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S, N) model.E, model.nu = E, nu surface = make_surface(N) rms = tm.Statistics2D.computeSpectralRMSSlope(surface) surface /= rms solver = tm.PolonskyKeerRey(model, surface, 1e-15) solver.solve(load) seq_tractions = np.copy(model.traction) seq_area = tm.Statistics2D.contact(model.traction) rms = comm.bcast(rms, root=0) model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S, N) model.E, model.nu = E, nu surface = make_surface(N) / rms solver = tm.PolonskyKeerRey(model, surface, 1e-15) solver.solve(load) tractions = model['traction'] area = tm.Statistics2D.contact(tractions) recv = comm.gather(tractions, root=0) if comm.rank == 0: tractions = np.concatenate(recv) error = np.linalg.norm(seq_tractions - tractions) / seq_tractions.size if False: print(error) import matplotlib.pyplot as plt plt.imshow(seq_tractions - tractions) plt.colorbar() plt.show() assert error < 1e-13 assert np.abs(seq_area - area) / seq_area < 1e-13 def model_for_dump(): model_type = tm.model_type.volume_2d discretization = [10, 2, 10] flat_domain = [1, 1] system_size = [0.5] + flat_domain # Creation of model model = tm.ModelFactory.createModel(model_type, system_size, discretization) model['displacement'][:] = MPI.COMM_WORLD.rank model['traction'][:] = MPI.COMM_WORLD.rank return model def mpi_h5dump(): try: from tamaas.dumpers import H5Dumper as Dumper import h5py except ImportError: return if not h5py.get_config().mpi: print('Did not test H5Dumper with MPI') return os.chdir(MPI.COMM_WORLD.bcast(tempfile.mkdtemp(), root=0)) print(os.getcwd()) model = model_for_dump() dumper_helper = Dumper('test_mpi_dump', 'displacement', 'traction') dumper_helper << model with h5py.File('hdf5/test_mpi_dump_0000.h5', 'r', driver='mpio', comm=MPI.COMM_WORLD) as handle: assert np.all(handle['displacement'][:, 0, :] == 0) assert np.all(handle['displacement'][:, 1, :] == 1) assert np.all(handle['traction'][0, :] == 0) assert np.all(handle['traction'][1, :] == 1) rmodel = dumper_helper.read('hdf5/test_mpi_dump_0000.h5') assert np.all(rmodel.displacement == MPI.COMM_WORLD.rank) assert np.all(rmodel.traction == MPI.COMM_WORLD.rank) def mpi_netcdfdump(): try: from tamaas.dumpers import NetCDFDumper as Dumper - import netCDF4 except ImportError: return - import subprocess - capture = subprocess.run(["nc-config", "--has-parallel4"], - capture_output=True) + try: + capture = subprocess.run(["nc-config", "--has-parallel4"], + capture_output=True) - if 'yes' not in capture.stdout.decode('ascii'): + if 'yes' not in capture.stdout.decode('ascii'): + print('Did not test NetCDF with MPI') + return + except FileNotFoundError: print('Did not test NetCDF with MPI') return os.chdir(MPI.COMM_WORLD.bcast(tempfile.mkdtemp(), root=0)) print(os.getcwd()) model = model_for_dump() dumper_helper = Dumper('test_mpi_dump', 'displacement', 'traction') dumper_helper << model rmodel = dumper_helper.read('netcdf/test_mpi_dump.nc') assert np.all(rmodel.displacement.astype(int) == MPI.COMM_WORLD.rank) assert np.all(rmodel.traction.astype(int) == MPI.COMM_WORLD.rank) def mpi_plastic_solve(): from conftest import UniformPlasticity from tamaas.nonlinear_solvers import DFSANECXXSolver as Solver patch_isotropic_plasticity = UniformPlasticity(tm.model_type.volume_2d, [1.] * 3, [4] * 3) model = patch_isotropic_plasticity.model residual = patch_isotropic_plasticity.residual applied_pressure = 0.1 solver = Solver(residual) solver.tolerance = 1e-15 pressure = model['traction'][..., 2] pressure[:] = applied_pressure solver.solve() solver.updateState() solution, normal = patch_isotropic_plasticity.solution(applied_pressure) for key in solution: error = norm(model[key] - solution[key]) / normal[key] assert error < 2e-15 def mpi_flood_fill(): contact = np.zeros([10, 1], dtype=np.bool) contact[:1, :] = True contact[-1:, :] = True clusters = tm.FloodFill.getClusters(contact, False) if tm.mpi.rank() == 0: assert len(clusters) == 3, "Wrong number of clusters" assert [c.getArea() for c in clusters] == [1, 2, 1], "Wrong areas" if __name__ == '__main__': mpi_polonsky_keer_compare()