Page MenuHomec4science
Files F62283274

mpi_routines.py
No OneTemporary
Actions

File Metadata

Created
Sun, May 12, 04:10

mpi_routines.py
View Options

# -*- coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-2020 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 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()
generator.setSizes(N)
generator.setSpectrum(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.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()
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.getDiscretization(), model.getBoundaryDiscretization()
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_size(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
if __name__ == '__main__':
mpi_polonsky_keer()

Event Timeline

c4science · Help