221 lines
10 KiB
Python
221 lines
10 KiB
Python
"""Sparse scientific frames and the approved new-reproduce presentation idiom.
|
|
|
|
Frame gears/UI colors/rotation sectors derive from DynamisLab continuous_video.py
|
|
(read-only SHA in docs). Modern velocity decoding and moving coordinates are used.
|
|
All camera limits, vorticity limits, and motion scale are fixed, never following y.
|
|
"""
|
|
from functools import lru_cache
|
|
import json
|
|
from pathlib import Path
|
|
import subprocess
|
|
import time
|
|
|
|
import numpy as np
|
|
|
|
D,U,CY=40.,.005,510.5
|
|
SIZE=(1280,800)
|
|
CW,CCW='#8854B8','#138F78'
|
|
MAIN=(0,120,1280,440)
|
|
INSET=(960,492,1240,772)
|
|
MAIN_LIMITS=(6.4,57.6,-6.4,6.4)
|
|
INSET_LIMITS=(28.75,33.75,-2.5,2.5)
|
|
|
|
|
|
def frame_schedule(steps,base_interval,accelerated=False):
|
|
"""Original 15-frame gear holds; only readback cadence changes."""
|
|
if accelerated:
|
|
ramp=15*(1+2+5)*base_interval
|
|
if steps<=2*ramp:
|
|
raise ValueError('too short for both display ramps')
|
|
blocks=[(s,15*s*base_interval) for s in (1,2,5)]
|
|
blocks += [(10,steps-2*ramp)]
|
|
blocks += [(s,15*s*base_interval) for s in (5,2,1)]
|
|
else:
|
|
blocks=[(1,steps)]
|
|
result,done=[],0
|
|
for speed,duration in blocks:
|
|
end=done+duration
|
|
while done<end:
|
|
done=min(done+speed*base_interval,end)
|
|
result.append((done,speed))
|
|
return result
|
|
|
|
|
|
def video_schedule(phase_steps,smoke=False):
|
|
if smoke:
|
|
phase=frame_schedule(phase_steps,8)
|
|
else:
|
|
# First 80 D/U of each phase is fast display; last 80 is normal.
|
|
half=phase_steps//2
|
|
phase=frame_schedule(half,800,True)
|
|
phase += [(half+n,s) for n,s in frame_schedule(phase_steps-half,800)]
|
|
return {p*phase_steps+n:s for p in range(3) for n,s in phase}
|
|
|
|
|
|
def vorticity(ux,uy,flags):
|
|
"""Centered dx=1 curl; mask every stencil touching nonfluid, no filling."""
|
|
ux,uy=np.asarray(ux,dtype=float),np.asarray(uy,dtype=float)
|
|
fluid=(flags&1)!=0
|
|
valid=(fluid[1:-1,1:-1]&fluid[1:-1,:-2]&fluid[1:-1,2:]&
|
|
fluid[:-2,1:-1]&fluid[2:,1:-1])
|
|
curl=((uy[1:-1,2:]-uy[1:-1,:-2])-(ux[2:,1:-1]-ux[:-2,1:-1]))/2
|
|
result=np.full(ux.shape,np.nan)
|
|
result[1:-1,1:-1]=np.where(valid,curl*D/U,np.nan)
|
|
return result
|
|
|
|
|
|
def pixel(x,y,box,limits):
|
|
left,top,right,bottom=box
|
|
x0,x1,y0,y1=limits
|
|
return left+(x-x0)/(x1-x0)*(right-left),bottom-(y-y0)/(y1-y0)*(bottom-top)
|
|
|
|
|
|
@lru_cache(maxsize=12)
|
|
def font(size,bold=False):
|
|
from PIL import ImageFont
|
|
name='DejaVuSans-Bold.ttf' if bold else 'DejaVuSans.ttf'
|
|
return ImageFont.truetype('/usr/share/fonts/truetype/dejavu/'+name,size)
|
|
|
|
|
|
def draw_circle(draw,x,y,r,wall_speed,box,limits):
|
|
cx,cy=pixel(x,y,box,limits)
|
|
radius=r*(box[2]-box[0])/(limits[1]-limits[0])
|
|
rect=(cx-radius,cy-radius,cx+radius,cy+radius)
|
|
draw.ellipse(rect,fill='white',outline='#263D4B',width=2)
|
|
if wall_speed is None:
|
|
return
|
|
fraction=abs(wall_speed)/(12*U)
|
|
if not np.isfinite(fraction) or fraction>1.00001:
|
|
raise RuntimeError('applied rotation exceeds fixed display reference')
|
|
r2=radius*.72
|
|
inner=(cx-r2,cy-r2,cx+r2,cy+r2)
|
|
draw.ellipse(inner,fill='#E4E9ED')
|
|
angle=min(fraction,1)*360
|
|
start,end=(-90,-90+angle) if wall_speed<0 else (-90-angle,-90)
|
|
if angle>0:
|
|
draw.pieslice(inner,start,end,fill=CW if wall_speed<0 else CCW)
|
|
|
|
|
|
class Video:
|
|
def __init__(self,out,scene,bodies,phase_steps,smoke=False,source='unknown'):
|
|
from PIL import Image
|
|
from matplotlib import colormaps
|
|
self.out,self.scene,self.bodies=out,scene,bodies
|
|
self.schedule=video_schedule(phase_steps,smoke)
|
|
self.phase_steps=phase_steps
|
|
self.source=source
|
|
self.records=[]
|
|
self.palette=np.rint(colormaps['RdBu_r'](np.linspace(0,1,21))[:,:3]*255).astype('uint8')
|
|
self.log=open(out/'ffmpeg.log','wb')
|
|
self.proc=subprocess.Popen(['ffmpeg','-hide_banner','-loglevel','warning','-n',
|
|
'-f','rawvideo','-pix_fmt','rgb24','-s','1280x800','-r','30','-i','pipe:0',
|
|
'-an','-c:v','libx264','-threads','4','-preset','fast','-crf','18',
|
|
'-pix_fmt','yuv420p','-movflags','+faststart',str(out/'continuous_viv.mp4')],
|
|
stdin=subprocess.PIPE,stderr=self.log)
|
|
self.closed=False
|
|
self.seconds=0.
|
|
|
|
def capture(self,sim,step,phase,eff,rows):
|
|
speed=self.schedule.get(step)
|
|
if speed is None:
|
|
return
|
|
from PIL import Image,ImageDraw
|
|
start=time.monotonic()
|
|
macro=sim.get_macroscopic()
|
|
flags=sim.get_flags().reshape(sim.field.ny,sim.field.nx)
|
|
fluid=(flags&1)!=0
|
|
if any(not np.isfinite(macro[k][fluid]).all() for k in ('rho','ux','uy')):
|
|
raise RuntimeError('nonfinite fluid field at video sample')
|
|
omega=vorticity(macro['ux'],macro['uy'],flags)
|
|
finite=np.isfinite(omega)
|
|
index=np.floor((np.clip(np.nan_to_num(omega),-3,3)+3)/6*21).astype(int).clip(0,20)
|
|
rgb=self.palette[index]
|
|
rgb[~finite]=[232,235,237]
|
|
field=Image.fromarray(rgb[::-1])
|
|
canvas=Image.new('RGB',SIZE,'white')
|
|
draw=ImageDraw.Draw(canvas)
|
|
for box,limits in ((MAIN,MAIN_LIMITS),(INSET,INSET_LIMITS)):
|
|
x0,x1,y0,y1=limits
|
|
# Lattice cell-centre coordinates -> image pixel-edge coordinates.
|
|
crop=(x0*D+.5,sim.field.ny-(CY+y1*D+.5),
|
|
x1*D+.5,sim.field.ny-(CY+y0*D+.5))
|
|
picture=field.crop(crop).resize((box[2]-box[0],box[3]-box[1]),Image.Resampling.BOX)
|
|
canvas.paste(picture,box[:2])
|
|
if box==INSET:
|
|
_,eq=pixel(30,0,box,limits)
|
|
for x in range(box[0],box[2],10):
|
|
draw.line((x,eq,min(x+5,box[2]),eq),fill='#657887',width=1)
|
|
if self.scene=='karman':
|
|
draw_circle(draw,10,0,1,None,box,limits) if box==MAIN else None
|
|
for bid,w in zip(self.bodies,eff):
|
|
obj=sim.bodies.objects[bid]
|
|
x,y=obj.motion.centre(obj.motion.applied)
|
|
draw_circle(draw,x/D,(y-CY)/D,obj.radius/D,float(w),box,limits)
|
|
x0,y0=pixel(INSET_LIMITS[0],INSET_LIMITS[2],MAIN,MAIN_LIMITS)
|
|
x1,y1=pixel(INSET_LIMITS[1],INSET_LIMITS[3],MAIN,MAIN_LIMITS)
|
|
draw.rectangle((x0,y1,x1,y0),outline='#647C8C',width=1)
|
|
draw.rectangle(INSET,outline='#647C8C',width=1)
|
|
def text(x,y,value,size=16,color='#263D4B',bold=False):
|
|
draw.text((x,y),value,font=font(size,bold),fill=color)
|
|
on=phase=='on'
|
|
text(28,15,'INCOMING FLOW',12,'#667887',True)
|
|
text(28,39,'Steady inflow' if self.scene=='steady' else 'Periodic cylinder wake',29,bold=True)
|
|
text(28,83,'One continuous flow and oscillator trajectory | Re_D = 50 | zero structural damping',15)
|
|
text(850,18,'CYLINDER ROTATION',12,bold=True)
|
|
text(850,43,'Clockwise',14,CW);text(1020,43,'Counterclockwise',14,CCW)
|
|
text(850,70,'Sector angle = actual surface speed',13,'#667887')
|
|
text(850,90,'Full circle = 12 x reference inflow speed',13,'#667887')
|
|
draw.line((24,453,1256,453),fill='#E5EBEF')
|
|
text(28,465,'CLOAKING CONTROL',12,'#667887',True)
|
|
draw.rounded_rectangle((28,490,92,524),radius=17,fill='#137F73' if on else '#A8B4BC')
|
|
cx=75 if on else 45
|
|
draw.ellipse((cx-12,495,cx+12,519),fill='white')
|
|
text(108,491,'ON' if on else 'OFF',22,'#137F73' if on else '#667887',True)
|
|
text(195,493,('Fixed rotations' if self.scene=='steady' else 'Frozen PPO') if on else 'Zero rotation',18)
|
|
text(570,467,f't U/D = {step*U/D:.2f}',18,bold=True)
|
|
local=(step-1)%self.phase_steps+1
|
|
transition=(phase!='off_before' and local*U/D<=1.)
|
|
text(570,495,f'Playback {speed}x | '+('switch transient' if transition else phase.replace('_',' ')),16)
|
|
text(960,467,f'Fixed lab-frame inset | y/D={rows[-1,1]/D:+.3f}',13,bold=True)
|
|
# Motion trace: fixed +/-D limits, a fixed 60 D/U trailing time span.
|
|
box=(55,560,895,732)
|
|
draw.rectangle(box,outline='#CCD6DD')
|
|
middle=(box[1]+box[3])/2
|
|
draw.line((box[0],middle,box[2],middle),fill='#A8B4BC')
|
|
t=step*U/D
|
|
left=max(0.,t-60)
|
|
first=max(0,int(left*D/U))
|
|
sample=rows[first::max(1,(len(rows)-first)//1200)]
|
|
xy=[(box[0]+(r[0]*U/D-left)/60*(box[2]-box[0]),middle-r[1]/D*(box[3]-box[1])/2) for r in sample]
|
|
if len(xy)>1:draw.line(xy,fill='#263D4B',width=2)
|
|
text(28,535,'Transverse displacement y/D (fixed scale; dashed inset line = equilibrium)',14)
|
|
text(24,554,'+1',12);text(25,642,'0',12);text(24,722,'-1',12)
|
|
text(55,739,f'{left:.1f}',12);text(835,739,f'{left+60:.1f}',12)
|
|
text(350,739,'response time t U/D',12)
|
|
text(28,774,f'Vorticity D/U: fixed [-3, +3] | exploratory, not convergence proof | source {self.source[:8]}',12)
|
|
self.proc.stdin.write(canvas.tobytes())
|
|
local=step%self.phase_steps
|
|
if not self.records or local in (0,self.phase_steps//2) or self.records[-1]['phase']!=phase:
|
|
canvas.save(self.out/f'keyframe-{step:09d}.png')
|
|
elapsed=time.monotonic()-start
|
|
self.seconds+=elapsed
|
|
self.records.append(dict(frame=len(self.records),step=step,absolute_clock=sim.stepper.step_count,
|
|
tU_D=t,phase=phase,playback_speed=speed,y_D=float(rows[-1,1]/D),
|
|
surface_speeds=eff.tolist(),render_s=elapsed,
|
|
fraction_finite_outside_color_limits=float(np.mean(np.abs(omega[finite])>3))))
|
|
|
|
def close(self):
|
|
if self.closed:return
|
|
self.closed=True
|
|
self.proc.stdin.close()
|
|
code=self.proc.wait(timeout=60)
|
|
self.log.close()
|
|
receipt=dict(frames=self.records,frame_count=len(self.records),fps=30,size=SIZE,
|
|
duration_s=len(self.records)/30,render_total_s=self.seconds,ffmpeg_exit=code,
|
|
normalization='centered dx=1 vorticity D/U, all 5 stencil cells fluid',
|
|
color_limits=[-3,3],color_bins=21,main_limits_D=MAIN_LIMITS,inset_limits_D=INSET_LIMITS,
|
|
y_origin_lattice=CY,trace_limits_D=[-1,1],inset_camera='fixed lab frame',
|
|
acceleration='frame cadence only; native solver and SI unchanged')
|
|
(self.out/'video.json').write_text(json.dumps(receipt,indent=2)+'\n')
|
|
if code:raise RuntimeError(f'ffmpeg exit {code}')
|