Fr0zka/VoxelForge
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

experiment: LivingPlane readable + alive mode (L163)

Browse Source 
RGBA32F displacement/velocity/imprint field; rolling base relief, grazing light, elevation colours, 7.5 cm contours,
follow camera; scenes: roll, rest and wobble, lift-off, second sphere; alive mode = travelling swells + viscoelastic
imprint. Intel ADL GT2: ~0.17 ms/step at 512^2. Open: run-to-run difference with queries on/off.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016dhxM43dnWvUu2B5Ph1iG1
This commit is contained in:
Fr0zkaandClaude Opus 5.5 committed 2026-11-02 10:55:32 +01:00
1 parent 5b15595e4a
commit e1d87eed26
3 files changed
+460 -128
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

Experiments/LivingPlane/LivingPlane.cpp
+388 -127
View File
@@ -8,6 +8,7 @@
#include <algorithm>
#include <array>
#include <cerrno>
#include <chrono>
#include <cctype>
#include <cmath>
#include <cstdint>
@@ -34,18 +35,33 @@ namespace
constexpr double DampingPerSecond=1.35;
constexpr double FoundationStiffnessPerSecondSquared=28.0;
constexpr double SphereRadiusMeters=0.43;
constexpr double SphereSigmaMeters=0.46;
constexpr double SpherePushMetersPerSecondSquared=14.0;
constexpr double SphereSigmaMeters=0.58;
constexpr double SpherePushMetersPerSecondSquared=36.0;
constexpr GLuint TimerQueryRingSize=64;
constexpr double Pi=3.14159265358979323846;
struct FBodyState
{
double X=0.0,Y=0.0,Z=0.0,Force=0.0;
};
struct FOptions
{
std::string Device="auto";
std::filesystem::path Output="Saved/Linux/l162/living-plane";
std::filesystem::path Output="Saved/Linux/l163/living-plane";
std::string Scene="wake";
int Size=128;
int Steps=360;
int CaptureEvery=1;
bool TimerQueries=true;
bool Breathing=false;
bool BreathingConfigured=false;
bool Imprint=false;
bool ImprintConfigured=false;
double BreathingStrength=4.8;
double ImprintCaptureMetersPerSecond=0.15;
double ImprintRecoverySeconds=6.0;
double ImprintMaxDepthMeters=0.34;
};
struct FDeviceCandidate
@@ -361,8 +377,8 @@ namespace
{
static constexpr char Source[]=R"GLSL(#version 430 core
layout(local_size_x=16,local_size_y=16) in;
layout(binding=0,rg32f) readonly uniform image2D uSource;
layout(binding=1,rg32f) writeonly uniform image2D uTarget;
layout(binding=0,rgba32f) readonly uniform image2D uSource;
layout(binding=1,rgba32f) writeonly uniform image2D uTarget;
uniform int uGridSize;
uniform float uDt;
uniform float uDx;
@@ -371,38 +387,61 @@ uniform float uWaveSpeed;
uniform float uDamping;
uniform float uFoundationStiffness;
uniform float uTime;
uniform vec2 uForcePosition;
uniform vec2 uForcePositions[2];
uniform float uForceSigma;
uniform float uForceStrength;
uniform float uForceStrengths[2];
uniform int uBreathing;
uniform float uBreathingStrength;
uniform int uImprint;
uniform float uImprintCaptureRate;
uniform float uImprintRecoverySeconds;
uniform float uImprintMaxDepth;
void main()
{
ivec2 p=ivec2(gl_GlobalInvocationID.xy);
if(p.x>=uGridSize||p.y>=uGridSize)return;
ivec2 lo=ivec2(0);
ivec2 hi=ivec2(uGridSize-1);
float h=imageLoad(uSource,p).r;
float v=imageLoad(uSource,p).g;
vec4 state=imageLoad(uSource,p);
float h=state.r;
float v=state.g;
float imprint=state.b;
float l=imageLoad(uSource,ivec2(max(p.x-1,lo.x),p.y)).r;
float r=imageLoad(uSource,ivec2(min(p.x+1,hi.x),p.y)).r;
float d=imageLoad(uSource,ivec2(p.x,max(p.y-1,lo.y))).r;
float u=imageLoad(uSource,ivec2(p.x,min(p.y+1,hi.y))).r;
float laplacian=(l+r+d+u-4.0*h)/(uDx*uDx);
vec2 world=(vec2(p)+vec2(0.5))*uDx-vec2(uPlaneLength*0.5);
vec2 offset=world-uForcePosition;
float sigma2=max(uForceSigma*uForceSigma,1.0e-8);
float interaction=-uForceStrength*exp(-dot(offset,offset)/(2.0*sigma2));
float interaction=0.0;
float contact=0.0;
for(int i=0;i<2;++i)
{
vec2 offset=world-uForcePositions[i];
float weight=exp(-dot(offset,offset)/(2.0*sigma2));
interaction-=uForceStrengths[i]*weight;
if(uForceStrengths[i]>0.0)contact=max(contact,weight);
}
float breathing=0.0;
if(uBreathing!=0)
{
float spatial=sin(world.x*0.72)*sin(world.y*0.72);
breathing=2.4*spatial*sin(uTime*1.45);
float swell=sin(world.x*0.48-uTime*0.22)
*cos(world.y*0.42+uTime*0.16);
breathing=uBreathingStrength*swell*sin(uTime*0.48);
}
float acceleration=uWaveSpeed*uWaveSpeed*laplacian-uFoundationStiffness*h
float nextImprint=imprint;
if(uImprint!=0)
{
nextImprint=clamp(imprint*exp(-uDt/max(uImprintRecoverySeconds,1.0e-4))
-contact*uImprintCaptureRate*uDt,-uImprintMaxDepth,0.0);
}
float foundationTarget=uImprint!=0?imprint:0.0;
float acceleration=uWaveSpeed*uWaveSpeed*laplacian
-uFoundationStiffness*(h-foundationTarget)
+interaction+breathing;
float nextVelocity=(v+acceleration*uDt)*exp(-uDamping*uDt);
float nextHeight=h+nextVelocity*uDt;
imageStore(uTarget,p,vec4(nextHeight,nextVelocity,0.0,0.0));
imageStore(uTarget,p,vec4(nextHeight,nextVelocity,nextImprint,0.0));
})GLSL";
const GLuint Shader=CompileShader(GL_COMPUTE_SHADER,Source);
if(!Shader)return 0;
@@ -424,11 +463,17 @@ uniform float uPlaneLength;
uniform int uRenderSphere;
uniform vec3 uSphereCenter;
uniform float uSphereRadius;
uniform int uSphereMaterial;
uniform mat4 uViewProjection;
out vec3 vWorld;
out vec3 vNormal;
out float vHeight;
flat out int vSphere;
flat out int vSphereMaterial;
float TerrainHeight(vec2 p)
{
return 0.105*sin(p.x*0.62+p.y*0.13)*cos(p.y*0.57)
+0.045*sin(p.x*1.04-p.y*0.78);
}
void main()
{
vec3 world;
@@ -442,11 +487,21 @@ void main()
ivec2 right=ivec2(min(texel.x+1,uGridSize-1),texel.y);
ivec2 down=ivec2(texel.x,max(texel.y-1,0));
ivec2 up=ivec2(texel.x,min(texel.y+1,uGridSize-1));
height=texelFetch(uField,texel,0).r;
float sampleDx=uPlaneLength/float(uGridSize);
vec2 leftWorld=(vec2(left)+vec2(0.5))*sampleDx-vec2(uPlaneLength*0.5);
vec2 rightWorld=(vec2(right)+vec2(0.5))*sampleDx-vec2(uPlaneLength*0.5);
vec2 downWorld=(vec2(down)+vec2(0.5))*sampleDx-vec2(uPlaneLength*0.5);
vec2 upWorld=(vec2(up)+vec2(0.5))*sampleDx-vec2(uPlaneLength*0.5);
float leftHeight=texelFetch(uField,left,0).r+TerrainHeight(leftWorld);
float rightHeight=texelFetch(uField,right,0).r+TerrainHeight(rightWorld);
float downHeight=texelFetch(uField,down,0).r+TerrainHeight(downWorld);
float upHeight=texelFetch(uField,up,0).r+TerrainHeight(upWorld);
height=texelFetch(uField,texel,0).r+TerrainHeight(
(vec2(texel)+vec2(0.5))*sampleDx-vec2(uPlaneLength*0.5));
float dx=max(float(right.x-left.x)*uPlaneLength/float(uGridSize),1.0e-6);
float dy=max(float(up.y-down.y)*uPlaneLength/float(uGridSize),1.0e-6);
float dhdx=(texelFetch(uField,right,0).r-texelFetch(uField,left,0).r)/dx;
float dhdy=(texelFetch(uField,up,0).r-texelFetch(uField,down,0).r)/dy;
float dhdx=(rightHeight-leftHeight)/dx;
float dhdy=(upHeight-downHeight)/dy;
world=vec3((aUv-vec2(0.5))*uPlaneLength,height);
normal=normalize(vec3(-dhdx,-dhdy,1.0));
}
@@ -459,37 +514,43 @@ void main()
vWorld=world;
vNormal=normal;
vHeight=height;
vSphere=uRenderSphere;
vSphereMaterial=uRenderSphere==0?0:uSphereMaterial;
gl_Position=uViewProjection*vec4(world,1.0);
})GLSL";
static constexpr char FragmentSource[]=R"GLSL(#version 430 core
in vec3 vWorld;
in vec3 vNormal;
in float vHeight;
flat in int vSphere;
flat in int vSphereMaterial;
uniform vec3 uCameraPosition;
out vec4 outColor;
void main()
{
vec3 n=normalize(vNormal);
vec3 l=normalize(vec3(-0.45,-0.55,0.72));
vec3 l=normalize(vec3(-0.62,-0.54,0.32));
vec3 v=normalize(uCameraPosition-vWorld);
float diffuse=max(dot(n,l),0.0);
float spec=pow(max(dot(reflect(-l,n),v),0.0),vSphere==0?58.0:38.0);
if(vSphere!=0)
float slope=1.0-clamp(n.z,0.0,1.0);
float spec=pow(max(dot(reflect(-l,n),v),0.0),vSphereMaterial==0?72.0:38.0);
if(vSphereMaterial!=0)
{
vec3 base=vec3(0.98,0.43,0.13);
vec3 base=vSphereMaterial==1?vec3(0.96,0.34,0.10):vec3(0.18,0.70,0.92);
float light=0.24+0.84*diffuse;
outColor=vec4(base*light+vec3(0.65,0.43,0.22)*spec,1.0);
}
else
{
float crest=clamp(0.50+vHeight*0.8,0.0,1.0);
vec3 deep=vec3(0.035,0.22,0.34);
vec3 crestColor=vec3(0.13,0.62,0.70);
vec3 base=mix(deep,crestColor,crest);
float light=0.28+0.88*diffuse;
outColor=vec4(base*light+vec3(0.31,0.62,0.68)*spec,1.0);
float low= smoothstep(-0.46,-0.08,vHeight);
float high=smoothstep(0.02,0.32,vHeight);
vec3 base=mix(vec3(0.025,0.075,0.13),vec3(0.055,0.30,0.31),low);
base=mix(base,vec3(0.28,0.46,0.28),high*0.82);
float contourPhase=fract((vHeight+0.60)/0.075);
float contourDistance=min(contourPhase,1.0-contourPhase);
float contour=smoothstep(0.025,0.16,contourDistance);
base=mix(base*0.35,base,contour);
base+=vec3(0.20,0.22,0.12)*smoothstep(0.08,0.36,slope);
float light=0.34+0.98*diffuse;
outColor=vec4(base*light+vec3(0.30,0.52,0.48)*spec,1.0);
}
})GLSL";
const GLuint VS=CompileShader(GL_VERTEX_SHADER,VertexSource);
@@ -509,6 +570,65 @@ void main()
struct FVec3 { double X=0.0,Y=0.0,Z=0.0; };
double TerrainHeight(double X,double Y)
{
return 0.105*std::sin(X*0.62+Y*0.13)*std::cos(Y*0.57)
+0.045*std::sin(X*1.04-Y*0.78);
}
double Smooth01(double T)
{
T=std::clamp(T,0.0,1.0);
return T*T*(3.0-2.0*T);
}
std::array<FBodyState,2> GetBodiesAt(double Time)
{
constexpr double RollEnd=1.62;
constexpr double RestEnd=2.55;
constexpr double LiftEnd=2.95;
constexpr double SecondaryStart=2.70;
constexpr double SecondaryEnd=5.55;
FBodyState Primary;
Primary.Y=-0.15;
if(Time<RollEnd)
{
Primary.X=-3.55+4.05*Smooth01(Time/RollEnd);
Primary.Z=TerrainHeight(Primary.X,Primary.Y)+0.41;
Primary.Force=SpherePushMetersPerSecondSquared;
}
else if(Time<RestEnd)
{
const double RestTime=Time-RollEnd;
Primary.X=0.50;
Primary.Z=TerrainHeight(Primary.X,Primary.Y)+0.41
+0.065*std::sin(RestTime*2.0*Pi*1.15);
Primary.Force=SpherePushMetersPerSecondSquared
*(0.82+0.18*std::sin(RestTime*2.0*Pi*1.15));
}
else if(Time<LiftEnd)
{
const double Lift=Smooth01((Time-RestEnd)/(LiftEnd-RestEnd));
Primary.X=0.50;
Primary.Z=TerrainHeight(Primary.X,Primary.Y)+0.41+1.10*Lift;
}
else
{
Primary.X=0.50;
Primary.Z=TerrainHeight(Primary.X,Primary.Y)+1.51;
}
FBodyState Secondary;
Secondary.Y=1.05+0.14*std::sin(Time*0.72);
const double Travel=Smooth01((Time-SecondaryStart)/(SecondaryEnd-SecondaryStart));
Secondary.X=5.15-10.0*Travel;
Secondary.Z=TerrainHeight(Secondary.X,Secondary.Y)+0.41;
if(Time>=SecondaryStart&&Time<=SecondaryEnd)
Secondary.Force=SpherePushMetersPerSecondSquared*0.78;
return {Primary,Secondary};
}
FVec3 Sub(FVec3 A,FVec3 B){return {A.X-B.X,A.Y-B.Y,A.Z-B.Z};}
FVec3 Cross(FVec3 A,FVec3 B)
{return {A.Y*B.Z-A.Z*B.Y,A.Z*B.X-A.X*B.Z,A.X*B.Y-A.Y*B.X};}
@@ -519,9 +639,8 @@ void main()
return Length>1.0e-12?FVec3{V.X/Length,V.Y/Length,V.Z/Length}:FVec3{};
}
std::array<float,16> MakeViewProjection(int Width,int Height,FVec3 Camera)
std::array<float,16> MakeViewProjection(int Width,int Height,FVec3 Camera,FVec3 Target)
{
const FVec3 Target={0.0,0.0,0.0};
const FVec3 Forward=Normalize(Sub(Target,Camera));
const FVec3 Side=Normalize(Cross(Forward,{0.0,0.0,1.0}));
const FVec3 Up=Cross(Side,Forward);
@@ -534,8 +653,8 @@ void main()
View[10]=static_cast<float>(-Forward.Z);View[14]=static_cast<float>(Dot(Forward,Camera));
View[15]=1.0f;
const double Near=0.1,Far=100.0;
const double F=1.0/std::tan(42.0*Pi/360.0);
const double Near=0.1,Far=60.0;
const double F=1.0/std::tan(39.0*Pi/360.0);
const double Aspect=static_cast<double>(Width)/static_cast<double>(Height);
std::array<float,16> Projection{};
Projection[0]=static_cast<float>(F/Aspect);
@@ -635,7 +754,7 @@ void main()
GLsizei GridIndexCount=0,SphereIndexCount=0;
int Width=0,Height=0,GridSize=0;
std::array<float,16> ViewProjection{};
FVec3 Camera{8.8,-12.2,9.4};
FVec3 Camera{};
void Destroy()
{
@@ -672,7 +791,7 @@ void main()
if(glCheckFramebufferStatus(GL_FRAMEBUFFER)!=GL_FRAMEBUFFER_COMPLETE)
{std::cerr<<"[LivingPlane] offscreen framebuffer is incomplete\n";return false;}
const int Segments=std::min(Size,256);
const int Segments=std::min(Size,512);
std::vector<FGridVertex> GridVertices;
std::vector<std::uint32_t> GridIndices;
MakeGrid(Segments,GridVertices,GridIndices);
@@ -708,14 +827,16 @@ void main()
static_cast<GLsizeiptr>(SphereIndices.size()*sizeof(std::uint32_t)),
SphereIndices.data(),GL_STATIC_DRAW);
ViewProjection=MakeViewProjection(Width,Height,Camera);
glBindVertexArray(0);
return glGetError()==GL_NO_ERROR;
}
bool Render(GLuint FieldTexture,const std::filesystem::path& Output,
double SphereX,double SphereY,double SphereZ)
const FBodyState& Primary,const FBodyState& Secondary)
{
const FVec3 Target={Primary.X,Primary.Y,0.0};
Camera={Primary.X+2.65,Primary.Y-4.45,3.65};
ViewProjection=MakeViewProjection(Width,Height,Camera,Target);
glBindFramebuffer(GL_FRAMEBUFFER,Framebuffer);
glViewport(0,0,Width,Height);
glDisable(GL_CULL_FACE);
@@ -739,11 +860,18 @@ void main()
glBindVertexArray(GridVao);
glDrawElements(GL_TRIANGLES,GridIndexCount,GL_UNSIGNED_INT,nullptr);
glUniform1i(glGetUniformLocation(Program,"uRenderSphere"),1);
glUniform3f(glGetUniformLocation(Program,"uSphereCenter"),
static_cast<float>(SphereX),static_cast<float>(SphereY),static_cast<float>(SphereZ));
glUniform1f(glGetUniformLocation(Program,"uSphereRadius"),
static_cast<float>(SphereRadiusMeters));
glBindVertexArray(SphereVao);
glUniform1i(glGetUniformLocation(Program,"uSphereMaterial"),1);
glUniform3f(glGetUniformLocation(Program,"uSphereCenter"),
static_cast<float>(Primary.X),static_cast<float>(Primary.Y),
static_cast<float>(Primary.Z));
glDrawElements(GL_TRIANGLES,SphereIndexCount,GL_UNSIGNED_INT,nullptr);
glUniform1i(glGetUniformLocation(Program,"uSphereMaterial"),2);
glUniform3f(glGetUniformLocation(Program,"uSphereCenter"),
static_cast<float>(Secondary.X),static_cast<float>(Secondary.Y),
static_cast<float>(Secondary.Z));
glDrawElements(GL_TRIANGLES,SphereIndexCount,GL_UNSIGNED_INT,nullptr);
glFinish();
std::vector<std::uint8_t> Pixels(static_cast<std::size_t>(Width)*Height*4);
@@ -768,6 +896,18 @@ void main()
return true;
}
bool ParseDouble(std::string_view Text,double& Value)
{
if(Text.empty())return false;
char* End=nullptr;
errno=0;
const std::string Copy(Text);
const double Parsed=std::strtod(Copy.c_str(),&End);
if(errno||!End||*End!='\0'||!std::isfinite(Parsed))return false;
Value=Parsed;
return true;
}
bool ParseOptions(int Argc,char** Argv,FOptions& Options)
{
for(int I=1;I<Argc;++I)
@@ -782,13 +922,17 @@ void main()
std::string_view Value;
if(Arg=="--help"||Arg=="-h")
{
std::cout<<"Usage: living-plane [--device=auto|gpu|cpu] [--size=N] [--steps=N] "
"[--out=DIR] [--capture-every=N] [--breathing=on|off]\n"
" --size is the square height-field and PNG resolution.\n"
" --capture-every=0 writes only timing.json. Default is 1.\n";
std::cout<<"Usage: living-plane [--device=auto|gpu|cpu] [--scene=wake|alive] "
"[--size=N] [--steps=N] [--out=DIR] [--capture-every=N]\n"
" [--queries=on|off] [--breathing=on|off] [--imprint=on|off]\n"
" [--breathing-strength=MPS2] [--imprint-capture=MPS] "
"[--imprint-recovery=SECONDS] [--imprint-depth=M]\n"
" --capture-every=0 writes only timing.json. --queries=off measures "
"batched dispatch wall time without per-step query overhead.\n";
std::exit(0);
}
else if(Take("--device",Value))Options.Device=Value;
else if(Take("--scene",Value))Options.Scene=Value;
else if(Take("--size",Value))
{if(!ParseInt(Value,Options.Size)){std::cerr<<"invalid --size\n";return false;}}
else if(Take("--steps",Value))
@@ -796,22 +940,66 @@ void main()
else if(Take("--capture-every",Value))
{if(!ParseInt(Value,Options.CaptureEvery)){std::cerr<<"invalid --capture-every\n";return false;}}
else if(Take("--out",Value))Options.Output=Value;
else if(Take("--queries",Value))
{
if(Value=="on")Options.TimerQueries=true;
else if(Value=="off")Options.TimerQueries=false;
else {std::cerr<<"--queries must be on or off\n";return false;}
}
else if(Take("--breathing",Value))
{
if(Value=="on")Options.Breathing=true;
else if(Value=="off")Options.Breathing=false;
else {std::cerr<<"--breathing must be on or off\n";return false;}
Options.BreathingConfigured=true;
}
else if(Take("--imprint",Value))
{
if(Value=="on")Options.Imprint=true;
else if(Value=="off")Options.Imprint=false;
else {std::cerr<<"--imprint must be on or off\n";return false;}
Options.ImprintConfigured=true;
}
else if(Take("--breathing-strength",Value))
{
if(!ParseDouble(Value,Options.BreathingStrength))
{std::cerr<<"invalid --breathing-strength\n";return false;}
}
else if(Take("--imprint-capture",Value))
{
if(!ParseDouble(Value,Options.ImprintCaptureMetersPerSecond))
{std::cerr<<"invalid --imprint-capture\n";return false;}
}
else if(Take("--imprint-recovery",Value))
{
if(!ParseDouble(Value,Options.ImprintRecoverySeconds))
{std::cerr<<"invalid --imprint-recovery\n";return false;}
}
else if(Take("--imprint-depth",Value))
{
if(!ParseDouble(Value,Options.ImprintMaxDepthMeters))
{std::cerr<<"invalid --imprint-depth\n";return false;}
}
else {std::cerr<<"unknown option: "<<Arg<<"\n";return false;}
}
if(Options.Device!="auto"&&Options.Device!="gpu"&&Options.Device!="cpu")
{std::cerr<<"--device must be auto, gpu, or cpu\n";return false;}
if(Options.Scene!="wake"&&Options.Scene!="alive")
{std::cerr<<"--scene must be wake or alive\n";return false;}
if(!Options.BreathingConfigured)Options.Breathing=Options.Scene=="alive";
if(!Options.ImprintConfigured)Options.Imprint=Options.Scene=="alive";
if(Options.Size<16||Options.Size>4096)
{std::cerr<<"--size must be in [16,4096]\n";return false;}
if(Options.Steps<1||Options.Steps>100000)
{std::cerr<<"--steps must be in [1,100000]\n";return false;}
if(Options.CaptureEvery<0||Options.CaptureEvery>100000)
{std::cerr<<"--capture-every must be in [0,100000]\n";return false;}
if(Options.BreathingStrength<0.0||Options.BreathingStrength>30.0
||Options.ImprintCaptureMetersPerSecond<0.0
||Options.ImprintCaptureMetersPerSecond>1.0
||Options.ImprintRecoverySeconds<0.1||Options.ImprintRecoverySeconds>60.0
||Options.ImprintMaxDepthMeters<0.0||Options.ImprintMaxDepthMeters>1.0)
{std::cerr<<"alive-mode parameters are outside their supported ranges\n";return false;}
return true;
}
@@ -826,8 +1014,9 @@ void main()
bool WriteTimingJson(const FOptions& Options,const FContext& Context,
std::string_view Renderer,std::string_view Vendor,std::string_view Version,
double Dt,const std::vector<GLuint64>& Times,const std::vector<double>& Heights,
const std::vector<double>& Velocities,std::size_t Frames,bool FirstStepAudit,
bool DynamicsAudit)
const std::vector<double>& Velocities,std::size_t Frames,int TimerQueryBits,
int TimerQueryWaits,bool TimerResultsAvailable,double DispatchWallMsPerStep,
bool FirstStepAudit,bool DynamicsAudit)
{
std::error_code Error;
std::filesystem::create_directories(Options.Output,Error);
@@ -862,6 +1051,7 @@ void main()
<<" \"gl_version\": \""<<JsonEscape(Version)<<"\",\n"
<<" \"size\": "<<Options.Size<<",\n"
<<" \"steps\": "<<Options.Steps<<",\n"
<<" \"scene\": \""<<Options.Scene<<"\",\n"
<<" \"plane_length_m\": "<<PlaneLengthMeters<<",\n"
<<" \"dt_seconds\": "<<Dt<<",\n"
<<" \"spectral_step_ratio\": "<<SpectralStep<<",\n"
@@ -869,12 +1059,30 @@ void main()
<<" \"damping_per_s\": "<<DampingPerSecond<<",\n"
<<" \"foundation_stiffness_per_s2\": "<<FoundationStiffnessPerSecondSquared<<",\n"
<<" \"breathing\": "<<(Options.Breathing?"true":"false")<<",\n"
<<" \"breathing_strength_m_per_s2\": "<<Options.BreathingStrength<<",\n"
<<" \"imprint\": "<<(Options.Imprint?"true":"false")<<",\n"
<<" \"imprint_capture_m_per_s\": "<<Options.ImprintCaptureMetersPerSecond<<",\n"
<<" \"imprint_recovery_seconds\": "<<Options.ImprintRecoverySeconds<<",\n"
<<" \"imprint_max_depth_m\": "<<Options.ImprintMaxDepthMeters<<",\n"
<<" \"frames_written\": "<<Frames<<",\n"
<<" \"first_step_equation_audit\": \""<<(FirstStepAudit?"pass":"fail")<<"\",\n"
<<" \"field_dynamics_audit\": \""<<(DynamicsAudit?"pass":"fail")<<"\",\n"
<<" \"timing_interpretation\": \""
<<(Context.ExecutionDevice=="gpu"?"GPU GL_TIME_ELAPSED per compute dispatch"
:"software-renderer GL_TIME_ELAPSED values; not GPU performance timings")<<"\",\n"
<<" \"timer_queries_enabled\": "<<(Options.TimerQueries?"true":"false")<<",\n"
<<" \"timer_query_ring_size\": "<<(Options.TimerQueries
?std::min<std::size_t>(TimerQueryRingSize,Times.size()):0)<<",\n"
<<" \"timer_query_counter_bits\": "<<TimerQueryBits<<",\n"
<<" \"timer_query_results_available_after_completion\": ";
if(Options.TimerQueries)Json<<(TimerResultsAvailable?"true":"false");
else Json<<"null";
Json<<",\n"
<<" \"timer_query_waits_before_reuse\": "<<TimerQueryWaits<<",\n"
<<" \"dispatch_wall_ms_per_step\": ";
if(DispatchWallMsPerStep>=0.0)Json<<DispatchWallMsPerStep;
else Json<<"null";
Json<<",\n \"dispatch_wall_window\": \""
<<(Options.CaptureEvery==0
?"simulation loop plus final glFinish; no frame capture; first-step equation readback included"
:"not suitable for step cost because frame capture is enabled")<<"\",\n"
<<" \"field_min_m\": "<<HeightMin<<",\n"
<<" \"field_max_m\": "<<HeightMax<<",\n"
<<" \"field_mean_m\": "<<HeightMean<<",\n"
@@ -924,7 +1132,7 @@ void main()
WaveSpeedMetersPerSecond/(Dx*Dx)+FoundationStiffnessPerSecondSquared);
std::array<GLuint,2> Fields{};
glGenTextures(2,Fields.data());
const std::vector<float> Zeroes(static_cast<std::size_t>(Size)*Size*2,0.0f);
const std::vector<float> Zeroes(static_cast<std::size_t>(Size)*Size*4,0.0f);
for(const GLuint Field:Fields)
{
glBindTexture(GL_TEXTURE_2D,Field);
@@ -932,17 +1140,20 @@ void main()
glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER,GL_NEAREST);
glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_WRAP_S,GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_WRAP_T,GL_CLAMP_TO_EDGE);
glTexImage2D(GL_TEXTURE_2D,0,GL_RG32F,Size,Size,0,GL_RG,GL_FLOAT,Zeroes.data());
glTexImage2D(GL_TEXTURE_2D,0,GL_RGBA32F,Size,Size,0,GL_RGBA,GL_FLOAT,Zeroes.data());
}
const int CaptureWidth=Size,CaptureHeight=Size;
FRenderResources Render;
if(Options.CaptureEvery>0&&!Render.Initialize(Size,CaptureWidth,CaptureHeight))
if(Options.CaptureEvery>0&&!Render.Initialize(Size,Size,Size))
{Render.Destroy();glDeleteTextures(2,Fields.data());glDeleteProgram(ComputeProgram);return false;}
const GLuint QueryCount=static_cast<GLuint>(Options.Steps);
std::vector<GLuint> Queries(QueryCount);
glGenQueries(static_cast<GLsizei>(QueryCount),Queries.data());
const int QueryCount=Options.TimerQueries
?std::min<int>(static_cast<int>(TimerQueryRingSize),Options.Steps):0;
std::vector<GLuint> Queries(static_cast<std::size_t>(QueryCount));
if(QueryCount>0)glGenQueries(QueryCount,Queries.data());
GLint TimerQueryBits=0;
if(Options.TimerQueries)glGetQueryiv(GL_TIME_ELAPSED,GL_QUERY_COUNTER_BITS,&TimerQueryBits);
const GLint LocationGrid=glGetUniformLocation(ComputeProgram,"uGridSize");
const GLint LocationDt=glGetUniformLocation(ComputeProgram,"uDt");
const GLint LocationDx=glGetUniformLocation(ComputeProgram,"uDx");
@@ -951,14 +1162,20 @@ void main()
const GLint LocationDamping=glGetUniformLocation(ComputeProgram,"uDamping");
const GLint LocationFoundation=glGetUniformLocation(ComputeProgram,"uFoundationStiffness");
const GLint LocationTime=glGetUniformLocation(ComputeProgram,"uTime");
const GLint LocationForcePosition=glGetUniformLocation(ComputeProgram,"uForcePosition");
const GLint LocationForcePositions=glGetUniformLocation(ComputeProgram,"uForcePositions[0]");
const GLint LocationForceSigma=glGetUniformLocation(ComputeProgram,"uForceSigma");
const GLint LocationForceStrength=glGetUniformLocation(ComputeProgram,"uForceStrength");
const GLint LocationForceStrengths=glGetUniformLocation(ComputeProgram,"uForceStrengths[0]");
const GLint LocationBreathing=glGetUniformLocation(ComputeProgram,"uBreathing");
const GLint LocationBreathingStrength=glGetUniformLocation(ComputeProgram,"uBreathingStrength");
const GLint LocationImprint=glGetUniformLocation(ComputeProgram,"uImprint");
const GLint LocationImprintCapture=glGetUniformLocation(ComputeProgram,"uImprintCaptureRate");
const GLint LocationImprintRecovery=glGetUniformLocation(ComputeProgram,"uImprintRecoverySeconds");
const GLint LocationImprintDepth=glGetUniformLocation(ComputeProgram,"uImprintMaxDepth");
if(LocationGrid<0||LocationDt<0||LocationDx<0||LocationPlane<0||LocationWaveSpeed<0
||LocationDamping<0||LocationFoundation<0||LocationTime<0
||LocationForcePosition<0||LocationForceSigma<0
||LocationForceStrength<0||LocationBreathing<0)
||LocationForcePositions<0||LocationForceSigma<0||LocationForceStrengths<0
||LocationBreathing<0||LocationBreathingStrength<0||LocationImprint<0
||LocationImprintCapture<0||LocationImprintRecovery<0||LocationImprintDepth<0)
{std::cerr<<"[LivingPlane] missing compute shader uniforms\n";return false;}
glUseProgram(ComputeProgram);
@@ -971,19 +1188,30 @@ void main()
glUniform1f(LocationFoundation,static_cast<float>(FoundationStiffnessPerSecondSquared));
glUniform1f(LocationForceSigma,static_cast<float>(SphereSigmaMeters));
glUniform1i(LocationBreathing,Options.Breathing?1:0);
glUniform1f(LocationBreathingStrength,static_cast<float>(Options.BreathingStrength));
glUniform1i(LocationImprint,Options.Imprint?1:0);
glUniform1f(LocationImprintCapture,
static_cast<float>(Options.ImprintCaptureMetersPerSecond));
glUniform1f(LocationImprintRecovery,static_cast<float>(Options.ImprintRecoverySeconds));
glUniform1f(LocationImprintDepth,static_cast<float>(Options.ImprintMaxDepthMeters));
int Current=0;
const auto Dispatch=[&](int Target,double Time,double SphereX,double SphereY,
double ForceStrength,GLuint Query)
const auto Dispatch=[&](int Target,double Time,
const std::array<FBodyState,2>& Bodies,GLuint Query)
{
glUseProgram(ComputeProgram);
glUniform1f(LocationTime,static_cast<float>(Time));
glUniform2f(LocationForcePosition,static_cast<float>(SphereX),static_cast<float>(SphereY));
glUniform1f(LocationForceStrength,static_cast<float>(ForceStrength));
const std::array<float,4> Positions={
static_cast<float>(Bodies[0].X),static_cast<float>(Bodies[0].Y),
static_cast<float>(Bodies[1].X),static_cast<float>(Bodies[1].Y)};
const std::array<float,2> Strengths={
static_cast<float>(Bodies[0].Force),static_cast<float>(Bodies[1].Force)};
glUniform2fv(LocationForcePositions,2,Positions.data());
glUniform1fv(LocationForceStrengths,2,Strengths.data());
glBindImageTexture(0,Fields[static_cast<std::size_t>(Current)],0,GL_FALSE,0,
GL_READ_ONLY,GL_RG32F);
GL_READ_ONLY,GL_RGBA32F);
glBindImageTexture(1,Fields[static_cast<std::size_t>(Target)],0,GL_FALSE,0,
GL_WRITE_ONLY,GL_RG32F);
GL_WRITE_ONLY,GL_RGBA32F);
if(Query)glBeginQuery(GL_TIME_ELAPSED,Query);
glDispatchCompute(static_cast<GLuint>((Size+15)/16),
static_cast<GLuint>((Size+15)/16),1);
@@ -992,14 +1220,14 @@ void main()
Current=Target;
};
// Warm the runtime-compiled compute kernel with zero forcing before timing.
for(int Warm=0;Warm<8;++Warm)
Dispatch(1-Current,-8.0*Dt,0.0,0.0,0.0,0);
// Warm the runtime-compiled compute kernel before measurement.
const std::array<FBodyState,2> NoBodies{};
for(int Warm=0;Warm<8;++Warm)Dispatch(1-Current,-8.0*Dt,NoBodies,0);
glFinish();
for(const GLuint Field:Fields)
{
glBindTexture(GL_TEXTURE_2D,Field);
glTexSubImage2D(GL_TEXTURE_2D,0,0,0,Size,Size,GL_RG,GL_FLOAT,Zeroes.data());
glTexSubImage2D(GL_TEXTURE_2D,0,0,0,Size,Size,GL_RGBA,GL_FLOAT,Zeroes.data());
}
glMemoryBarrier(GL_TEXTURE_UPDATE_BARRIER_BIT|GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);
Current=0;
@@ -1007,62 +1235,72 @@ void main()
FRenderResources* RenderPtr=Options.CaptureEvery>0?&Render:nullptr;
std::size_t FramesWritten=0;
int FrameIndex=0;
auto SaveFrame=[&](double SphereX,double SphereY,double SphereZ)->bool
auto SaveFrame=[&](double Time)->bool
{
if(!RenderPtr)return true;
char Filename[32]{};
std::snprintf(Filename,sizeof(Filename),"frame_%06d.png",FrameIndex++);
const auto Bodies=GetBodiesAt(Time);
if(!Render.Render(Fields[static_cast<std::size_t>(Current)],
Options.Output/Filename,SphereX,SphereY,SphereZ))return false;
Options.Output/Filename,Bodies[0],Bodies[1]))return false;
++FramesWritten;
return true;
};
const int MoveSteps=std::max(2,static_cast<int>(std::lround(Options.Steps*0.58)));
const auto SphereAt=[&](int Step)
{
const double MoveT=std::clamp(static_cast<double>(Step)/std::max(1,MoveSteps-1),0.0,1.0);
const double Smooth=MoveT*MoveT*(3.0-2.0*MoveT);
const double X=-4.8+9.6*Smooth;
const bool Rest=Step>=MoveSteps;
const double RestTime=Rest?static_cast<double>(Step-MoveSteps)*Dt:0.0;
const double Bob=Rest?std::sin(RestTime*2.0*Pi*1.25):0.0;
const double Z=0.36+0.08*Bob;
const double Force=SpherePushMetersPerSecondSquared*(1.0-0.32*Bob);
return std::array<double,4>{X,0.0,Z,Force};
};
const auto Initial=SphereAt(0);
if(!SaveFrame(Initial[0],Initial[1],Initial[2]))
{Render.Destroy();glDeleteQueries(static_cast<GLsizei>(Queries.size()),Queries.data());
if(!SaveFrame(0.0))
{Render.Destroy();if(QueryCount>0)glDeleteQueries(QueryCount,Queries.data());
glDeleteTextures(2,Fields.data());glDeleteProgram(ComputeProgram);return false;}
std::vector<GLuint64> Times(Options.TimerQueries
?static_cast<std::size_t>(Options.Steps):0);
int TimerQueryWaits=0;
bool TimerResultsAvailable=Options.TimerQueries;
auto ReadTimerQuery=[&](int QueryStep)
{
const GLuint Query=Queries[static_cast<std::size_t>(QueryStep%QueryCount)];
GLint Available=GL_FALSE;
glGetQueryObjectiv(Query,GL_QUERY_RESULT_AVAILABLE,&Available);
if(Available!=GL_TRUE)++TimerQueryWaits;
glGetQueryObjectui64v(Query,GL_QUERY_RESULT,
&Times[static_cast<std::size_t>(QueryStep)]);
glGetQueryObjectiv(Query,GL_QUERY_RESULT_AVAILABLE,&Available);
TimerResultsAvailable=TimerResultsAvailable&&Available==GL_TRUE;
};
std::vector<float> FirstField;
bool FirstStepAudit=false;
const auto WallStart=std::chrono::steady_clock::now();
for(int Step=0;Step<Options.Steps;++Step)
{
const auto Sphere=SphereAt(Step);
if(QueryCount>0&&Step>=QueryCount)ReadTimerQuery(Step-QueryCount);
const double Time=static_cast<double>(Step)*Dt;
const auto Bodies=GetBodiesAt(Time);
const int Target=1-Current;
Dispatch(Target,static_cast<double>(Step)*Dt,Sphere[0],Sphere[1],Sphere[3],
Queries[static_cast<std::size_t>(Step)]);
const GLuint Query=QueryCount>0
?Queries[static_cast<std::size_t>(Step%QueryCount)]:0;
Dispatch(Target,Time,Bodies,Query);
if(Step==0)
{
glFinish();
FirstField.resize(static_cast<std::size_t>(Size)*Size*2);
FirstField.resize(static_cast<std::size_t>(Size)*Size*4);
glBindTexture(GL_TEXTURE_2D,Fields[static_cast<std::size_t>(Current)]);
glGetTexImage(GL_TEXTURE_2D,0,GL_RG,GL_FLOAT,FirstField.data());
glGetTexImage(GL_TEXTURE_2D,0,GL_RGBA,GL_FLOAT,FirstField.data());
const int X=std::clamp(static_cast<int>(std::floor(
(Sphere[0]+PlaneLengthMeters*0.5)/Dx)),0,Size-1);
(Bodies[0].X+PlaneLengthMeters*0.5)/Dx)),0,Size-1);
const int Y=std::clamp(static_cast<int>(std::floor(
(Sphere[1]+PlaneLengthMeters*0.5)/Dx)),0,Size-1);
(Bodies[0].Y+PlaneLengthMeters*0.5)/Dx)),0,Size-1);
const double CellX=(X+0.5)*Dx-PlaneLengthMeters*0.5;
const double CellY=(Y+0.5)*Dx-PlaneLengthMeters*0.5;
const double Distance2=(CellX-Sphere[0])*(CellX-Sphere[0])
+(CellY-Sphere[1])*(CellY-Sphere[1]);
const double Force=-Sphere[3]*std::exp(-Distance2/
(2.0*SphereSigmaMeters*SphereSigmaMeters));
double Force=0.0;
for(const FBodyState& Body:Bodies)
{
const double Distance2=(CellX-Body.X)*(CellX-Body.X)
+(CellY-Body.Y)*(CellY-Body.Y);
Force-=Body.Force*std::exp(-Distance2/
(2.0*SphereSigmaMeters*SphereSigmaMeters));
}
const double Expected=Force*Dt*Dt*std::exp(-DampingPerSecond*Dt);
const std::size_t Index=(static_cast<std::size_t>(Y)*Size+X)*2;
const std::size_t Index=(static_cast<std::size_t>(Y)*Size+X)*4;
const double Actual=FirstField[Index];
const double Tolerance=std::max(2.0e-7,std::abs(Expected)*2.5e-3);
FirstStepAudit=std::isfinite(Actual)&&std::abs(Actual-Expected)<=Tolerance;
@@ -1073,54 +1311,77 @@ void main()
if(RenderPtr&&(Options.CaptureEvery==1
||((Step+1)%Options.CaptureEvery==0)||Step==Options.Steps-1))
{
const auto FrameSphere=SphereAt(Step);
if(!SaveFrame(FrameSphere[0],FrameSphere[1],FrameSphere[2]))
{Render.Destroy();glDeleteQueries(static_cast<GLsizei>(Queries.size()),Queries.data());
if(!SaveFrame(static_cast<double>(Step+1)*Dt))
{Render.Destroy();if(QueryCount>0)glDeleteQueries(QueryCount,Queries.data());
glDeleteTextures(2,Fields.data());glDeleteProgram(ComputeProgram);return false;}
}
}
glFinish();
const auto WallEnd=std::chrono::steady_clock::now();
const double WallMs=std::chrono::duration<double,std::milli>(WallEnd-WallStart).count();
const double DispatchWallMsPerStep=Options.CaptureEvery==0
?WallMs/Options.Steps:-1.0;
std::vector<GLuint64> Times(static_cast<std::size_t>(Options.Steps));
for(int Step=0;Step<Options.Steps;++Step)
glGetQueryObjectui64v(Queries[static_cast<std::size_t>(Step)],GL_QUERY_RESULT,
&Times[static_cast<std::size_t>(Step)]);
std::vector<float> FinalField(static_cast<std::size_t>(Size)*Size*2);
if(QueryCount>0)
for(int Step=std::max(0,Options.Steps-QueryCount);Step<Options.Steps;++Step)
ReadTimerQuery(Step);
const bool QueryAudit=!Options.TimerQueries||TimerResultsAvailable;
std::vector<float> FinalField(static_cast<std::size_t>(Size)*Size*4);
glBindTexture(GL_TEXTURE_2D,Fields[static_cast<std::size_t>(Current)]);
glGetTexImage(GL_TEXTURE_2D,0,GL_RG,GL_FLOAT,FinalField.data());
glGetTexImage(GL_TEXTURE_2D,0,GL_RGBA,GL_FLOAT,FinalField.data());
std::vector<double> Heights,Velocities;
Heights.reserve(static_cast<std::size_t>(Size)*Size);
Velocities.reserve(static_cast<std::size_t>(Size)*Size);
bool Finite=true;
for(std::size_t I=0;I<FinalField.size();I+=2)
double ImprintMin=0.0;
for(std::size_t I=0;I<FinalField.size();I+=4)
{
Finite=Finite&&std::isfinite(FinalField[I])&&std::isfinite(FinalField[I+1]);
Finite=Finite&&std::isfinite(FinalField[I])&&std::isfinite(FinalField[I+1])
&&std::isfinite(FinalField[I+2]);
Heights.push_back(FinalField[I]);
Velocities.push_back(FinalField[I+1]);
ImprintMin=std::min(ImprintMin,static_cast<double>(FinalField[I+2]));
}
const auto [HeightMinIt,HeightMaxIt]=std::minmax_element(Heights.begin(),Heights.end());
const bool DynamicsAudit=Finite&&*HeightMinIt<-1.0e-7&&*HeightMaxIt-*HeightMinIt<2.0;
const bool Good=FirstStepAudit&&DynamicsAudit&&glGetError()==GL_NO_ERROR;
const bool DynamicsAudit=Finite&&*HeightMinIt<-1.0e-7&&*HeightMaxIt-*HeightMinIt<2.0
&&(!Options.Imprint||ImprintMin>=-Options.ImprintMaxDepthMeters-1.0e-5);
const bool Good=FirstStepAudit&&DynamicsAudit&&QueryAudit&&glGetError()==GL_NO_ERROR;
if(!DynamicsAudit)
std::cerr<<"[LivingPlane] field audit failed: finite="<<Finite<<" min="
<<*HeightMinIt<<" max="<<*HeightMaxIt<<"\n";
<<*HeightMinIt<<" max="<<*HeightMaxIt<<" imprint_min="<<ImprintMin<<"\n";
if(!QueryAudit)
std::cerr<<"[LivingPlane] timer-query result availability audit failed after glFinish\n";
if(!WriteTimingJson(Options,Context,Renderer,Vendor,Version,Dt,Times,Heights,Velocities,
FramesWritten,FirstStepAudit,DynamicsAudit))
{Render.Destroy();glDeleteQueries(static_cast<GLsizei>(Queries.size()),Queries.data());
FramesWritten,TimerQueryBits,TimerQueryWaits,TimerResultsAvailable,
DispatchWallMsPerStep,FirstStepAudit,DynamicsAudit))
{Render.Destroy();if(QueryCount>0)glDeleteQueries(QueryCount,Queries.data());
glDeleteTextures(2,Fields.data());glDeleteProgram(ComputeProgram);return false;}
double MeanNs=0.0;
for(const GLuint64 Ns:Times)MeanNs+=static_cast<double>(Ns);
MeanNs/=Times.size();
std::cout<<"[LivingPlane] size="<<Size<<" steps="<<Options.Steps<<" dt_ms="
<<Dt*1000.0<<" GL_TIME_ELAPSED_mean_ms="<<MeanNs*1.0e-6
<<" median_ms="<<Quantile([&]{std::vector<double> V;for(auto N:Times)V.push_back(N*1.0e-6);return V;}(),0.5)
<<" field_min_m="<<*HeightMinIt<<" field_max_m="<<*HeightMaxIt
<<" frames="<<FramesWritten<<" equation_audit="<<(FirstStepAudit?"PASS":"FAIL")
const std::vector<double> QueryMs=[&]
{
std::vector<double> Result;
Result.reserve(Times.size());
for(const GLuint64 Ns:Times)Result.push_back(static_cast<double>(Ns)*1.0e-6);
return Result;
}();
std::cout<<"[LivingPlane] scene="<<Options.Scene<<" size="<<Size
<<" steps="<<Options.Steps<<" dt_ms="<<Dt*1000.0
<<" dispatch_wall_ms_per_step=";
if(DispatchWallMsPerStep>=0.0)std::cout<<DispatchWallMsPerStep;
else std::cout<<"n/a(capture enabled)";
if(Options.TimerQueries)
std::cout<<" GL_TIME_ELAPSED_median_ms="<<Quantile(QueryMs,0.50)
<<" timer_query_bits="<<TimerQueryBits<<" timer_query_waits="<<TimerQueryWaits;
std::cout<<" field_min_m="<<*HeightMinIt<<" field_max_m="<<*HeightMaxIt
<<" imprint_min_m="<<ImprintMin<<" frames="<<FramesWritten
<<" equation_audit="<<(FirstStepAudit?"PASS":"FAIL")
<<" timer_query_audit="<<(Options.TimerQueries
?(QueryAudit?"PASS":"FAIL"):"NOT_RUN")
<<" status="<<(Good?"PASS":"FAIL")<<" out="<<Options.Output.string()<<"\n";
Render.Destroy();
glDeleteQueries(static_cast<GLsizei>(Queries.size()),Queries.data());
if(QueryCount>0)glDeleteQueries(QueryCount,Queries.data());
glDeleteTextures(2,Fields.data());
glDeleteProgram(ComputeProgram);
return Good;