Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/mpi_routines.py
9578f56e3c51features/netcdf-mpi

Last Change
Raw File

mpi_routines.py
View Options

# -*- coding: utf-8 -*-
# @file
# LICENSE
#
# Copyright (©) 2016-2021 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 <https://www.gnu.org/licenses/>.
from __future__ import print_function
import shutil
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 = 0
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_n, local_offset = tm.mpi.local_shape(N), tm.mpi.local_offset(N)
local_surface = np.copy(
surface[local_offset:local_offset+local_n[0], :])
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 = solve()
precv = comm.gather(tractions, root=0)
drecv = comm.gather(displacements, root=0)
if comm.rank == 0:
tractions = np.concatenate(precv)
displacements = np.concatenate(drecv)
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'])
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']
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
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
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
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)
if tm.mpi.rank() == 0:
shutil.rmtree('hdf5')
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
if __name__ == '__main__':
mpi_polonsky_keer_compare()
c4science · Help