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