Piston Shock#
import mmf_setup
mmf_setup.nbinit()
from gpe.imports import *
from importlib import reload
from gpe.Examples import piston
from gpe import minimize
This cell adds /home/docs/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:
- Choose "Trust Notebook" from the "File" menu.
- Re-execute this cell.
- Reload the notebook.
[I 20:01:29 numexpr.utils] NumExpr defaulting to 2 threads.
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[2], line 2
1 mmf_setup.nbinit()
----> 2 from gpe.imports import *
3 from importlib import reload
4 from gpe.Examples import piston
5 from gpe import minimize
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/imports.py:28
26 from mmfutils.plot import imcontourf # noqa: E402
27 from gpe.minimize import MinimizeState # noqa: E402
---> 28 from gpe.utils import evolve_to, evolve, evolves # noqa: E402
29 from gpe.plot_utils import MPLGrid # noqa: E402
30 from pytimeode.evolvers import EvolverSplit, EvolverABM # noqa: E402
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/utils.py:35
32 from persist.objects import Archivable
33 from persist.archive import Archive
---> 35 from pytimeode.evolvers import EvolverABM
36 from pytimeode.mixins import ArrayStateMixin
37 from pytimeode.interfaces import implementer, IStateForABMEvolvers
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/__init__.py:2
1 from . import interfaces
----> 2 from . import mixins
3 from . import evolvers
5 __all__ = ["interfaces", "mixins", "evolvers"]
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/mixins.py:39
30 from .interfaces import (
31 IState,
32 IStateApply,
(...) 36 implementer,
37 )
38 from . import interfaces
---> 39 from .utils import expr
41 __all__ = [
42 "StateMixin",
43 "StatesMixin",
(...) 47 "ArraysStateWithBraketMixin",
48 ]
51 @implementer(IState)
52 class StateMixin:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/utils/expr.py:18
15 import numpy as np
17 try:
---> 18 import sympy
19 except ImportError:
20 sympy = None
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/__init__.py:77
70 from .logic import (to_cnf, to_dnf, to_nnf, And, Or, Not, Xor, Nand, Nor,
71 Implies, Equivalent, ITE, POSform, SOPform, simplify_logic, bool_map,
72 true, false, satisfiable)
74 from .assumptions import (AppliedPredicate, Predicate, AssumptionsContext,
75 assuming, Q, ask, register_handler, remove_handler, refine)
---> 77 from .polys import (Poly, PurePoly, poly_from_expr, parallel_poly_from_expr,
78 degree, total_degree, degree_list, LC, LM, LT, pdiv, prem, pquo,
79 pexquo, div, rem, quo, exquo, half_gcdex, gcdex, invert,
80 subresultants, resultant, discriminant, cofactors, gcd_list, gcd,
81 lcm_list, lcm, terms_gcd, trunc, monic, content, primitive, compose,
82 decompose, sturm, gff_list, gff, sqf_norm, sqf_part, sqf_list, sqf,
83 factor_list, factor, intervals, refine_root, count_roots, all_roots,
84 real_roots, nroots, ground_roots, nth_power_roots_poly, cancel,
85 reduced, groebner, is_zero_dimensional, GroebnerBasis, poly,
86 symmetrize, horner, interpolate, rational_interpolate, viete, together,
87 BasePolynomialError, ExactQuotientFailed, PolynomialDivisionFailed,
88 OperationNotSupported, HeuristicGCDFailed, HomomorphismFailed,
89 IsomorphismFailed, ExtraneousFactors, EvaluationFailed,
90 RefinementFailed, CoercionFailed, NotInvertible, NotReversible,
91 NotAlgebraic, DomainError, PolynomialError, UnificationFailed,
92 GeneratorsError, GeneratorsNeeded, ComputationFailed,
93 UnivariatePolynomialError, MultivariatePolynomialError,
94 PolificationFailed, OptionError, FlagError, minpoly,
95 minimal_polynomial, primitive_element, field_isomorphism,
96 to_number_field, isolate, round_two, prime_decomp, prime_valuation,
97 galois_group, itermonomials, Monomial, lex, grlex,
98 grevlex, ilex, igrlex, igrevlex, CRootOf, rootof, RootOf,
99 ComplexRootOf, RootSum, roots, Domain, FiniteField, IntegerRing,
100 RationalField, RealField, ComplexField, PythonFiniteField,
101 GMPYFiniteField, PythonIntegerRing, GMPYIntegerRing, PythonRational,
102 GMPYRationalField, AlgebraicField, PolynomialRing, FractionField,
103 ExpressionDomain, FF_python, FF_gmpy, ZZ_python, ZZ_gmpy, QQ_python,
104 QQ_gmpy, GF, FF, ZZ, QQ, ZZ_I, QQ_I, RR, CC, EX, EXRAW,
105 construct_domain, swinnerton_dyer_poly, cyclotomic_poly,
106 symmetric_poly, random_poly, interpolating_poly, jacobi_poly,
107 chebyshevt_poly, chebyshevu_poly, hermite_poly, hermite_prob_poly,
108 legendre_poly, laguerre_poly, apart, apart_list, assemble_partfrac_list,
109 Options, ring, xring, vring, sring, field, xfield, vfield, sfield)
111 from .series import (Order, O, limit, Limit, gruntz, series, approximants,
112 residue, EmptySequence, SeqPer, SeqFormula, sequence, SeqAdd, SeqMul,
113 fourier_series, fps, difference_delta, limit_seq)
115 from .functions import (factorial, factorial2, rf, ff, binomial,
116 RisingFactorial, FallingFactorial, subfactorial, carmichael,
117 fibonacci, lucas, motzkin, tribonacci, harmonic, bernoulli, bell, euler,
(...) 138 Znm, elliptic_k, elliptic_f, elliptic_e, elliptic_pi, beta, mathieus,
139 mathieuc, mathieusprime, mathieucprime, riemann_xi, betainc, betainc_regularized)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/__init__.py:79
3 __all__ = [
4 'Poly', 'PurePoly', 'poly_from_expr', 'parallel_poly_from_expr', 'degree',
5 'total_degree', 'degree_list', 'LC', 'LM', 'LT', 'pdiv', 'prem', 'pquo',
(...) 65 'field', 'xfield', 'vfield', 'sfield'
66 ]
68 from .polytools import (Poly, PurePoly, poly_from_expr,
69 parallel_poly_from_expr, degree, total_degree, degree_list, LC, LM,
70 LT, pdiv, prem, pquo, pexquo, div, rem, quo, exquo, half_gcdex, gcdex,
(...) 76 nth_power_roots_poly, cancel, reduced, groebner, is_zero_dimensional,
77 GroebnerBasis, poly)
---> 79 from .polyfuncs import (symmetrize, horner, interpolate,
80 rational_interpolate, viete)
82 from .rationaltools import together
84 from .polyerrors import (BasePolynomialError, ExactQuotientFailed,
85 PolynomialDivisionFailed, OperationNotSupported, HeuristicGCDFailed,
86 HomomorphismFailed, IsomorphismFailed, ExtraneousFactors,
(...) 91 MultivariatePolynomialError, PolificationFailed, OptionError,
92 FlagError)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/polyfuncs.py:10
8 from sympy.polys.polyoptions import allowed_flags, build_options
9 from sympy.polys.polytools import poly_from_expr, Poly
---> 10 from sympy.polys.specialpolys import (
11 symmetric_poly, interpolating_poly)
12 from sympy.polys.rings import sring
13 from sympy.utilities import numbered_symbols, take, public
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/specialpolys.py:298
294 return dmp_mul(f, h, n, K), dmp_mul(g, h, n, K), h
296 # A few useful polynomials from Wang's paper ('78).
--> 298 from sympy.polys.rings import ring
300 def _f_0():
301 R, x, y, z = ring("x,y,z", ZZ)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/rings.py:31
27 from sympy.polys.polyoptions import (Domain as DomainOpt,
28 Order as OrderOpt, build_options)
29 from sympy.polys.polyutils import (expr_from_dict, _dict_reorder,
30 _parallel_dict_from_expr)
---> 31 from sympy.printing.defaults import DefaultPrinting
32 from sympy.utilities import public, subsets
33 from sympy.utilities.iterables import is_sequence
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/printing/__init__.py:11
7 from .mathml import mathml, print_mathml
9 from .python import python, print_python
---> 11 from .pycode import pycode
13 from .codeprinter import print_ccode, print_fcode
15 from .codeprinter import ccode, fcode, cxxcode, rust_code # noqa:F811
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/printing/pycode.py:11
9 from sympy.core.mod import Mod
10 from .precedence import precedence
---> 11 from .codeprinter import CodePrinter
13 _kw = {
14 'and', 'as', 'assert', 'break', 'class', 'continue', 'def', 'del', 'elif',
15 'else', 'except', 'finally', 'for', 'from', 'global', 'if', 'import', 'in',
16 'is', 'lambda', 'not', 'or', 'pass', 'raise', 'return', 'try', 'while',
17 'with', 'yield', 'None', 'False', 'nonlocal', 'True'
18 }
20 _known_functions = {
21 'Abs': 'abs',
22 'Min': 'min',
23 'Max': 'max',
24 }
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/printing/codeprinter.py:13
11 from sympy.core.sorting import default_sort_key
12 from sympy.core.symbol import Symbol
---> 13 from sympy.functions.elementary.complexes import re
14 from sympy.printing.str import StrPrinter
15 from sympy.printing.precedence import precedence, PRECEDENCE
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/functions/__init__.py:39
37 from sympy.functions.special.singularity_functions import SingularityFunction
38 from sympy.functions.special.delta_functions import DiracDelta, Heaviside
---> 39 from sympy.functions.special.bsplines import bspline_basis, bspline_basis_set, interpolating_spline
40 from sympy.functions.special.bessel import (besselj, bessely, besseli, besselk,
41 hankel1, hankel2, jn, yn, jn_zeros, hn1, hn2, airyai, airybi, airyaiprime, airybiprime, marcumq)
42 from sympy.functions.special.hyper import hyper, meijerg, appellf1
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/functions/special/bsplines.py:5
3 from sympy.functions import Piecewise, piecewise_fold
4 from sympy.logic.boolalg import And
----> 5 from sympy.sets.sets import Interval
7 from functools import lru_cache
10 def _ivl(cond, x):
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/sets/__init__.py:11
9 from .powerset import PowerSet
10 from ..core.singleton import S
---> 11 from .handlers.comparison import _eval_is_eq # noqa:F401
12 Complexes = S.Complexes
13 EmptySet = S.EmptySet
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/sets/handlers/comparison.py:46
42 eqs = (is_eq(x, y) for x, y in zip(lhs.sets, rhs.sets))
43 return tfn[fuzzy_and(map(fuzzy_bool, eqs))]
---> 46 @dispatch(Set, Basic) # type:ignore
47 def _eval_is_eq(lhs, rhs): # noqa: F811
48 return False
51 @dispatch(Set, Set) # type:ignore
52 def _eval_is_eq(lhs, rhs): # noqa: F811
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/multipledispatch/core.py:71, in dispatch.<locals>._(func)
68 namespace[name] = Dispatcher(name)
69 dispatcher = namespace[name]
---> 71 dispatcher.add(types, func, on_ambiguity=on_ambiguity)
72 return dispatcher
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/multipledispatch/dispatcher.py:219, in Dispatcher.add(self, signature, func, on_ambiguity)
213 raise TypeError("Tried to dispatch on non-type: %s\n"
214 "In signature: <%s>\n"
215 "In function: %s" %
216 (typ, str_sig, self.name))
218 self.funcs[signature] = func
--> 219 self.reorder(on_ambiguity=on_ambiguity)
220 self._cache.clear()
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/multipledispatch/dispatcher.py:225, in Dispatcher.reorder(self, on_ambiguity)
223 if _resolve[0]:
224 self.ordering = ordering(self.funcs)
--> 225 amb = ambiguities(self.funcs)
226 if amb:
227 on_ambiguity(self, amb)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/multipledispatch/conflict.py:27, in ambiguities(signatures)
25 """ All signature pairs such that A is ambiguous with B """
26 signatures = list(map(tuple, signatures))
---> 27 return {(a, b) for a in signatures for b in signatures
28 if hash(a) < hash(b)
29 and ambiguous(a, b)
30 and not any(supercedes(c, a) and supercedes(c, b)
31 for c in signatures)}
KeyboardInterrupt:
plt.rcParams["figure.figsize"] = (10, 5)
reload(piston)
e1 = piston.Experiment(tube=True)
# e2 = piston.Experiment(basis_type='axial')
# s = e.get_initial_state()
def callback(state, _n=[0]):
if _n[0] % 10 == 0:
plt.clf()
state.plot()
display(plt.gcf())
clear_output(wait=True)
_n[0] += 1
# s = minimize.MinimizeState(e1.get_state()).minimize(f_tol=0.1, callback=callback)
# s0 = e0.get_initial_state()
# s0.plot()
s1 = e1.get_initial_state()
# s2 = e2.get_state()
s1.plot()
s1.m
from mmfutils.plot import imcontourf
imcontourf(s2.xyz[0], s2.xyz[1], s2.get_density())
plt.xlim(-60, 250)
plt.ylim(-2, 2)
s.t = 0 * piston.u.ms
plt.plot(s.basis.xyz[0], s.experiment.get_Vt(s))
from mmfutils.plot import imcontourf
from gpe.plot_utils import MPLGrid
e = EvolverSplit(s1, dt=1.0 * s1.t_scale)
history = [s1.copy()]
e = EvolverSplit(history[-1], dt=0.5 * s1.t_scale)
with NoInterrupt() as interrupted:
while not interrupted and e.t < 120 * piston.u.ms:
e.evolve(1000)
history.append(e.get_y())
plt.clf()
# grid = MPLGrid()
s = e.y
# grid.grid()
s.plot()
# ax = grid.next()
# imcontourf(s.xyz[0], s.xyz[1], s.get_density())
# ax.set_xlim(0, 150)
# ax.set_aspect(1)
display(plt.gcf())
clear_output(wait=True)
data = np.asarray([_s[...] for _s in history])
np.save("piston.npy", data)
grid = MPLGrid()
s1.plot??
# (grid.grid())