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 03:38:05 numexpr.utils] NumExpr defaulting to 2 threads.
[I 03:38:05 root] Patching zope.interface.document.asReStructuredText to format code
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()
/home/docs/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/minimize.py:311: OptimizeWarning: Unknown solver options: tries
res = sp.optimize.minimize(
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[6], line 4
1 # s = minimize.MinimizeState(e1.get_state()).minimize(f_tol=0.1, callback=callback)
2 # s0 = e0.get_initial_state()
3 # s0.plot()
----> 4 s1 = e1.get_initial_state()
5 # s2 = e2.get_state()
6 s1.plot()
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/Examples/piston.py:237, in Experiment.get_initial_state(self, perturb, E_tol, psi_tol, disp, tries, cool_steps, cool_dt_t_scale, minimize, **kw)
235 if "use_scipy" not in kw:
236 kw["tries"] = tries
--> 237 state = m.minimize(E_tol=E_tol, psi_tol=psi_tol, disp=disp, **kw)
239 self._debug_state = state # Store in case evolve fails
241 if cool_steps > 1:
242 # Cool a bit to remove any fluctuations.
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/minimize.py:751, in MinimizeState.minimize(self, psi_tol, E_tol, callback, _debug, **kw)
749 return super().minimize(_debug=_debug, **kw)
750 else:
--> 751 x = super().minimize(**kw)
753 state = self.unpack(x)
755 if self.fix_N:
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/minimize.py:282, in Minimize.minimize(self, plot, callback, method, polish, broyden_alpha, broyden_opts, f_tol, x_tol, use_scipy, ignore_f, _test, _debug, _log, use_cache, bounds, **kw)
280 if Version(sp.__version__) > Version("1.15.0"):
281 options.pop("disp", None)
--> 282 res = self._minimize(
283 f=_f,
284 df=_df,
285 x0=_x[0],
286 method=method,
287 callback=callback_,
288 bounds=bounds,
289 options=options,
290 )
292 self.minimize_results = res
294 if not res.success:
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/minimize.py:311, in Minimize._minimize(self, f, df, x0, method, callback, bounds, options)
309 def _minimize(self, f, df, x0, method, callback, bounds, options):
310 """Interface to the scipy minimizer."""
--> 311 res = sp.optimize.minimize(
312 fun=f,
313 jac=df,
314 x0=x0,
315 method=method,
316 callback=callback,
317 bounds=bounds,
318 options=options,
319 )
320 res.f = f
321 res.df = df
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/optimize/_minimize.py:784, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
781 res = _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,
782 **options)
783 elif meth == 'l-bfgs-b':
--> 784 res = _minimize_lbfgsb(fun, x0, args, jac, bounds,
785 callback=callback, **options)
786 elif meth == 'tnc':
787 res = _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,
788 **options)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/optimize/_lbfgsb_py.py:469, in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, workers, **unknown_options)
461 _lbfgsb.setulb(m, x, low_bnd, upper_bnd, nbd, f, g, factr, pgtol, wa,
462 iwa, task, lsave, isave, dsave, maxls, ln_task)
464 if task[0] == 3:
465 # The minimization routine wants f and g at the current x.
466 # Note that interruptions due to maxfun are postponed
467 # until the completion of the current minimization iteration.
468 # Overwrite f and g:
--> 469 f, g = func_and_grad(x)
470 elif task[0] == 1:
471 # new iteration
472 n_iterations += 1
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:412, in ScalarFunction.fun_and_grad(self, x)
410 if not np.array_equal(x, self.x):
411 self._update_x(x)
--> 412 self._update_fun()
413 self._update_grad()
414 return self.f, self.g
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:362, in ScalarFunction._update_fun(self)
360 def _update_fun(self):
361 if not self.f_updated:
--> 362 fx = self._wrapped_fun(self.x)
363 self._nfev += 1
364 if fx < self._lowest_f:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/_lib/_util.py:603, in _ScalarFunctionWrapper.__call__(self, x)
600 def __call__(self, x):
601 # Send a copy because the user may overwrite it.
602 # The user of this class might want `x` to remain unchanged.
--> 603 fx = self.f(np.copy(x), *self.args)
604 self.nfev += 1
606 # Make sure the function returns a true scalar
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/minimize.py:158, in Minimize.minimize.<locals>._f(x)
152 if (
153 not use_cache
154 or _cache[1] is None
155 or not np.allclose(x, _cache[0], atol=1e-32, rtol=_EPS)
156 ):
157 _cache[0] = x.copy()
--> 158 _cache[1:] = self.f_df(x)
159 if _log:
160 self._calls.append(tuple(_cache))
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/minimize.py:663, in MinimizeState.f_df(self, x)
660 if self.fix_N:
661 s, N = psi.normalize()
--> 663 Hpsi = psi.get_Hy(subtract_mu=self.fix_N)
665 if self.fix_N:
666 Hpsi *= s
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/bec.py:877, in _StateBase.get_Hy(self, subtract_mu)
875 """Return `H(y)` for convenience only."""
876 dy = self.empty()
--> 877 self.compute_dy_dt(dy=dy, subtract_mu=subtract_mu)
878 Hy = dy / self._phase
879 return Hy
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/bec.py:775, in _StateBase.compute_dy_dt(self, dy, subtract_mu)
773 y = self
774 Ky = y.copy()
--> 775 Ky.apply_laplacian(factor=self.K_factor)
776 Vy = y.copy()
777 Vy.apply_V(V=self.get_V_GPU())
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/bec.py:1465, in StateScaleBase.apply_laplacian(self, factor, exp, **_kw)
1463 lams, dlams, ddlams = self.get_lambdas()
1464 k2 = sum((_k / _l) ** 2 for _l, _k in zip(lams, self.basis._pxyz))
-> 1465 super().apply_laplacian(factor=factor, exp=exp, k2=k2, **_kw)
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/bec.py:858, in _StateBase.apply_laplacian(self, factor, exp, **_kw)
856 if _v is not None:
857 _kw[_k] = _v
--> 858 self.data[...] = self.basis.laplacian(self.data, factor=factor, exp=exp, **_kw)
860 if self.Omega:
861 assert not exp
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/mmfutils/math/bases/bases.py:431, in PeriodicBasis.laplacian(self, y, factor, factors, exp, kx2, k2, kwz2, **_kw)
428 assert _raise_factors_err(factors=factors, kx2=kx2, k2=k2)
430 # Apply K
--> 431 yt = self.fftn(y)
432 laplacian_y = self.ifftn(
433 self._apply_K(yt, kx2=kx2, k2=k2, exp=exp, factor=factor, factors=factors)
434 )
436 if kwz2 != 0:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/mmfutils/math/bases/bases.py:500, in PeriodicBasis.fftn(self, x)
498 """Perform the fft along spatial axes"""
499 axes = self.axes % len(x.shape)
--> 500 return self._fftn(x, axes=axes)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/mmfutils/performance/fft.py:458, in fftn(a, s, axes)
456 if key not in _FFT_CACHE:
457 _FFT_CACHE[key] = get_fftn(a=a.copy(), s=s, axes=axes)
--> 458 res = _FFT_CACHE[key](a)
459 if _COPY_OUTPUT and not res.flags["OWNDATA"]:
460 res = res.copy()
KeyboardInterrupt:
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())