1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
#!/usr/bin/env python
# -----------------------------------------------------------------------------
# Programmer(s): David J. Gardner @ LLNL
# -----------------------------------------------------------------------------
# SUNDIALS Copyright Start
# Copyright (c) 2002-2022, Lawrence Livermore National Security
# and Southern Methodist University.
# All rights reserved.
#
# See the top-level LICENSE and NOTICE files for details.
#
# SPDX-License-Identifier: BSD-3-Clause
# SUNDIALS Copyright End
# -----------------------------------------------------------------------------
# matplotlib-based plotting script for cvsPraticle_dns.c example
# -----------------------------------------------------------------------------
# imports
import argparse
import numpy as np
import matplotlib.pyplot as plt
# command line options
parser = argparse.ArgumentParser(description='Plots cvsPraticle_dns output')
parser.add_argument('--sfile', type=str,
default='cvsParticle_solution.txt',
help='solution output file to read')
parser.add_argument('--efile', type=str,
default='cvsParticle_error.txt',
help='error output file to read')
parser.add_argument('--alpha', type=float, nargs=1,
default=1.0,
help='set a non-default alpha value')
parser.add_argument('--slim', type=float, nargs=2,
help='x and y limits for solution plot')
parser.add_argument('--eylim', type=float, nargs=2,
help='y limits for error plot')
# parse inputs
args = parser.parse_args()
# read solution output file
data = np.loadtxt(args.sfile, dtype=np.double)
# extract times and positions
t = data[:, 0]
x = data[:, 1]
y = data[:, 2]
# unit circle
tt = np.linspace(0,np.pi*2,10000)
xt = np.cos(tt)
yt = np.sin(tt)
# plot solution
fig, ax = plt.subplots()
plt.plot(xt, yt, color='black', linestyle='--')
plt.scatter(x, y, color='red')
if (args.slim):
plt.xlim((args.slim[0], args.slim[1]))
plt.ylim((args.slim[0], args.slim[1]))
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution')
ax.set_aspect('equal')
# true solution
xt = np.cos(args.alpha * t)
yt = np.sin(args.alpha * t)
# plot solution
fig, ax = plt.subplots()
plt.plot(t, x, linestyle='-', label='x')
plt.plot(t, xt, linestyle='--', label='x true')
plt.plot(t, y, linestyle='-', label='y')
plt.plot(t, yt, linestyle='--', label='y true')
plt.xlabel('t')
plt.ylabel('position')
plt.title('Particle Position Over Time')
plt.legend(loc='lower right')
# read error output file
data = np.loadtxt(args.efile, dtype=np.double)
# extract times, position errors, and constraint error
t = data[:, 0]
xerr = np.absolute(data[:, 1])
yerr = np.absolute(data[:, 2])
cerr = np.absolute(data[:, 3])
# plot solution
fig, ax = plt.subplots()
plt.semilogy(t, xerr, label='x err')
plt.semilogy(t, yerr, label='y err')
plt.semilogy(t, cerr, label='c err')
if (args.eylim):
plt.ylim((args.eylim[0], args.eylim[1]))
plt.xlabel('time')
plt.ylabel('error')
plt.legend(loc='lower right')
plt.title('Error in position and constraint')
plt.grid()
# display plots
plt.show()
##### end of script #####
|
