"""
Part of comet/pyhed/misc
"""
# COMET - COupled Magnetic resonance Electrical resistivity Tomography
# Copyright (C) 2019 Nico Skibbe
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU 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 General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
import numpy as np
import sys
import traceback
from contextlib import contextmanager
from comet.pyhed import log
try:
from mpi4py import MPI as mpi
except ImportError:
pass
[docs]@contextmanager
def abortIfError():
try:
yield
except:
try:
log.critical('MPI process {} tries to raise traceback.'
.format(mpi.COMM_WORLD.Get_rank()))
traceback.print_exc(file=sys.stdout)
sys.stdout.flush()
finally:
log.critical('MPI process {} about to exit due to critical error.'
.format(mpi.COMM_WORLD.Get_rank()))
mpi.COMM_WORLD.Abort(1)
[docs]def importCustemResults(name, ntx=1):
target_f = [df.Function(self.FS.V) for t in range(n_tx)]
#
f = df.HDF5File(self.MP.mpi_cs, f_name, "r")
f.read(target_f[ti], "data/vector_%d" % ti)
[docs]def saveFenicsField(savename_base, loop, secondary=False,
htr=None, hti=None, hsr=None, hsi=None):
"""
Save fenics fields from loppclass object in mpirun environment in
gimli single core sorting for later use in single core tasks.
savename total field :
- savename_base + '_total.npy'
if secondary is True:
savename secondary field:
- savename_base + '_secondary.npy'
Saved variables are:
- loop.secMOD.PP.H_t_r_cg,
- loop.secMOD.PP.H_t_i_cg,
- loop.secMOD.PP.H_s_r_cg,
- loop.secMOD.PP.H_s_i_cg
"""
from mpi4py import MPI
import pygimli as pg
from comet.pyhed.misc import convertCoordinates
import dolfin as df
com = MPI.COMM_WORLD
if loop.secMOD.FS.p == 2:
fs = df.VectorFunctionSpace(loop.secMOD.FS.mesh, 'CG', 1)
ht_real_cg = df.interpolate(loop.secMOD.PP.H_t_r_cg, fs)
hr_v = ht_real_cg.vector()
ht_imag_cg = df.interpolate(loop.secMOD.PP.H_t_i_cg, fs)
hi_v = ht_imag_cg.vector()
else:
# h total fields
if hasattr(loop.secMOD.PP, 'H_t_r_cg'):
if isinstance(loop.secMOD.PP.H_t_r_cg, list):
hr_v = loop.secMOD.PP.H_t_r_cg[0].vector()
hi_v = loop.secMOD.PP.H_t_i_cg[0].vector()
else:
hr_v = loop.secMOD.PP.H_t_r_cg.vector()
hi_v = loop.secMOD.PP.H_t_i_cg.vector()
else:
# no anomalies: export primary fields
if isinstance(loop.secMOD.PP.H_0_r_cg, list):
hr_v = loop.secMOD.PP.H_0_r_cg[0].vector()
hi_v = loop.secMOD.PP.H_0_i_cg[0].vector()
else:
hr_v = loop.secMOD.FS.PF.H_0_r_cg.vector()
hi_v = loop.secMOD.FS.PF.H_0_i_cg.vector()
# local distribution of global indices
r = hr_v.local_range()
# gather fields automatically
tot_real = hr_v.gather_on_zero().reshape(-1, 3)
tot_imag = hi_v.gather_on_zero().reshape(-1, 3)
if secondary is True:
if loop.secMOD.FS.p == 2:
if hasattr(loop.secMOD.PP, 'H_s_r'):
hs_real_cg = df.interpolate(loop.secMOD.PP.H_s_r, fs)
hr_vs = hs_real_cg.vector()
hs_imag_cg = df.interpolate(loop.secMOD.PP.H_s_i, fs)
hi_vs = hs_imag_cg.vector()
sec_real = hr_vs.gather_on_zero().reshape(-1, 3)
sec_imag = hi_vs.gather_on_zero().reshape(-1, 3)
else:
secondary = False
else:
if hasattr(loop.secMOD.PP, 'H_s_r_cg'):
if isinstance(loop.secMOD.PP.H_s_r_cg, list):
hr_vs = loop.secMOD.PP.H_s_r_cg[0].vector()
hi_vs = loop.secMOD.PP.H_s_i_cg[0].vector()
else:
hr_vs = loop.secMOD.PP.H_s_r_cg.vector()
hi_vs = loop.secMOD.PP.H_s_i_cg.vector()
sec_real = hr_vs.gather_on_zero().reshape(-1, 3)
sec_imag = hi_vs.gather_on_zero().reshape(-1, 3)
else:
secondary = False
# get distributed coords, fenics order
if loop.secMOD.FE.FS.p != 1:
V_cg = df.VectorFunctionSpace(loop.secMOD.FE.FS.mesh, 'CG', 1)
else:
V_cg = loop.secMOD.FE.FS.V_cg
coords_para = V_cg.tabulate_dof_coordinates().reshape(
-1, 3)[0::3, :]
# gather coords [ndarray, ndarray, ndarray], unsorted
coords_single = com.gather(coords_para, root=0)
# gather corresponding indices
ids_single = com.gather(r, root=0)
# %% continue only with one process
if com.Get_rank() == 0:
# gather coords in fenics, internally fenics sorted
df_c = np.zeros_like(tot_real).ravel()
for i, ids in enumerate(ids_single): # one run per process
df_c[ids[0]:ids[1]] = coords_single[i].ravel()
df_c = df_c.reshape(-1, 3)
# coords from gimli
gimliMesh = pg.Mesh('{}/_h5/{}.bms'.format(
loop.secondary_config.m_dir,
loop.secondary_config.mesh_name))
pg_c = gimliMesh.positions().array()
# calculation for resorting
log.debug('converted mesh')
pg_df, df_pg = convertCoordinates(pg_c, df_c)
# save numpy vec in gimli sorting, p1
totname = '{}_total.npy'.format(savename_base)
gimli_field_t = (tot_real + 1j * tot_imag)[df_pg].T
np.save(totname, gimli_field_t)
if secondary is True:
secname = '{}_secondary.npy'.format(savename_base)
gimli_field_s = (sec_real + 1j * sec_imag)[df_pg].T
np.save(secname, gimli_field_s)
log.info('Save secondary field: "{}"'.format(secname))
# optional vtk export
loop.loopmesh = gimliMesh
loop.field = gimli_field_t
loop.exportVTK('{}_total_pg_df.vtk'.format(savename_base))
if secondary is True:
loop.field = gimli_field_s
loop.exportVTK('{}_secondary_pg_df.vtk'.format(savename_base))
# The End
# to try stackoverflow magic
# !/usr/bin/env python
# from mpi4py import MPI
# import numpy
# import sys
# import os
# rank = MPI.COMM_WORLD.Get_rank()
# new_comm = MPI.COMM_WORLD.Split(color=rank, key=rank)
# print(new_comm.Get_rank())
# cwd=os.getcwd()
# os.mkdir(str(rank))
# directory=os.path.join(cwd,str(rank))
# print(rank,directory)
# os.chdir(directory)
#
#
# new_comm.Spawn("SOME_MPI_EXECUTABLE_HERE",
# args=[""],
# maxprocs=4)