Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/tests/mpi_routines.py
14f13ab72ad0maxwell_viscoelastic…

Last Change
Raw File

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
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 = [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.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-7
if __name__ == '__main__':
mpi_polonsky_keer()
c4science · Help