1#!/usr/bin/env python 2# ------------------------------------------------------------------------------ 3# Programmer(s): Daniel R. Reynolds @ SMU 4# David J. Gardner @ LLNL 5# ------------------------------------------------------------------------------ 6# SUNDIALS Copyright Start 7# Copyright (c) 2002-2021, Lawrence Livermore National Security 8# and Southern Methodist University. 9# All rights reserved. 10# 11# See the top-level LICENSE and NOTICE files for details. 12# 13# SPDX-License-Identifier: BSD-3-Clause 14# SUNDIALS Copyright End 15# ------------------------------------------------------------------------------ 16# matplotlib-based plotting script for the parallel cvheat2D_p examples 17# ------------------------------------------------------------------------------ 18 19# imports 20import sys, os 21import shlex 22import numpy as np 23from pylab import * 24from mpl_toolkits.mplot3d import Axes3D 25from matplotlib import cm 26import matplotlib.pyplot as plt 27 28# ------------------------------------------------------------------------------ 29 30# read MPI root process problem info file 31infofile = 'heat2d_info.00000.txt' 32 33with open(infofile) as fn: 34 35 # read the file line by line 36 for line in fn: 37 38 # split line into list 39 text = shlex.split(line) 40 41 # x-direction upper domian bound 42 if "xu" in line: 43 xu = float(text[1]) 44 continue 45 46 # y-direction upper domain bound 47 if "yu" in line: 48 yu = float(text[1]) 49 continue 50 51 # nodes in the x-direction (global) 52 if "nx" in line: 53 nx = int(text[1]) 54 continue 55 56 # nodes in the y-direction (global) 57 if "ny" in line: 58 ny = int(text[1]) 59 continue 60 61 # total number of MPI processes 62 if "np"in line: 63 nprocs = int(text[1]) 64 continue 65 66 # number of output times 67 if "nt" in line: 68 nt = int(text[1]) 69 continue 70 71# ------------------------------------------------------------------------------ 72 73# load subdomain information, store in table 74subdomains = np.zeros((nprocs,4), dtype=np.int) 75 76for i in range(nprocs): 77 78 infofile = 'heat2d_info.' + repr(i).zfill(5) + '.txt' 79 80 with open(infofile) as fn: 81 82 # read the file line by line 83 for line in fn: 84 85 # split line into list 86 text = shlex.split(line) 87 88 # x-direction starting index 89 if "is" in line: 90 subdomains[i,0] = float(text[1]) 91 continue 92 93 # x-direction ending index 94 if "ie" in line: 95 subdomains[i,1] = float(text[1]) 96 continue 97 98 # y-direction starting index 99 if "js" in line: 100 subdomains[i,2] = float(text[1]) 101 continue 102 103 # y-direction ending index 104 if "je" in line: 105 subdomains[i,3] = float(text[1]) 106 continue 107 108# ------------------------------------------------------------------------------ 109 110plottype = ['solution', 'error'] 111 112for pt in plottype: 113 114 # fill array with data 115 time = np.zeros(nt) 116 result = np.zeros((nt, ny, nx)) 117 118 for i in range(nprocs): 119 120 datafile = 'heat2d_' + pt + '.' + repr(i).zfill(5) + '.txt' 121 122 # load data 123 data = np.loadtxt(datafile, dtype=np.double) 124 125 if (np.shape(data)[0] != nt): 126 sys.exit('error: subdomain ' + i + ' has an incorrect number of time steps') 127 128 # subdomain indices 129 istart = subdomains[i,0] 130 iend = subdomains[i,1] 131 jstart = subdomains[i,2] 132 jend = subdomains[i,3] 133 nxl = iend - istart + 1 134 nyl = jend - jstart + 1 135 136 # extract data 137 for i in range(nt): 138 time[i] = data[i,0] 139 result[i,jstart:jend+1,istart:iend+1] = np.reshape(data[i,1:], (nyl,nxl)) 140 141 # determine extents of plots 142 maxtemp = 1.1 * result.max() 143 mintemp = 0.9 * result.min() 144 145 # set x and y meshgrid objects 146 xspan = np.linspace(0.0, xu, nx) 147 yspan = np.linspace(0.0, yu, ny) 148 X,Y = np.meshgrid(xspan, yspan) 149 150 nxstr = repr(nx) 151 nystr = repr(ny) 152 153 # generate plots 154 for tstep in range(nt): 155 156 # set string constants for output plots, current time, mesh size 157 pname = 'heat2d_surf_' + pt + '.' + repr(tstep).zfill(3) + '.png' 158 tstr = str(time[tstep]) 159 160 # plot surface and save to disk 161 fig = plt.figure(1) 162 ax = fig.add_subplot(111, projection='3d') 163 164 ax.plot_surface(X, Y, result[tstep,:,:], rstride=1, cstride=1, 165 cmap=cm.jet, linewidth=0, antialiased=True, shade=True) 166 167 ax.set_xlabel('x') 168 ax.set_ylabel('y') 169 ax.set_zlim((mintemp, maxtemp)) 170 ax.view_init(20,45) 171 if (pt == 'solution'): 172 title('u(x,y) at t = ' + tstr) 173 else: 174 title('error(x,y) at t = ' + tstr) 175 savefig(pname) 176 plt.close() 177 178##### end of script ##### 179