Geschwindichkeitsfeld divergenz frei machen
Verfasst: 30.07.2012, 18:52
Hallo zusammen,
ich habe jetzt schon seit Wochen das Problem das mein Fluid Solver nicht richtig funktioniert. Konkret benötigt man ja für die Massenerhaltung ein divergenzfreies Geschwindigkeitsfeld. Dazu Führe ich folgende Schritte im CS durch:
1.) Divergenz Vel berechnen:
2.) Jacob Iteration (x500):
3.) Projektion:
Zum testen vervende ich als Kraft nur die Gravitationskraft (0,-10,0) u. als Boundery einfach einen Würfel.
Wenn ich nach der Projektion das Div Feld neu berechnen ist leider nicht wirklich divergenzfrei.
Vielen Dank,
Manuel
ich habe jetzt schon seit Wochen das Problem das mein Fluid Solver nicht richtig funktioniert. Konkret benötigt man ja für die Massenerhaltung ein divergenzfreies Geschwindigkeitsfeld. Dazu Führe ich folgende Schritte im CS durch:
1.) Divergenz Vel berechnen:
Code: Alles auswählen
float3 get_quantity_bd(uint3 index)
{
float3 value;
if(index.x < 0)
return -q[index3_to_index(uint3(0,index.y,index.z))];
else if(index.x >= FLUID_GRID_CELL_COUNT.x)
return -q[index3_to_index(uint3(FLUID_GRID_CELL_COUNT.x-1,index.y,index.z))];
if(index.y < 0)
return -q[index3_to_index(uint3(index.x,0,index.z))];
else if(index.y >= FLUID_GRID_CELL_COUNT.y)
return -q[index3_to_index(uint3(index.x,FLUID_GRID_CELL_COUNT.y-1,index.z))];
if(index.z < 0)
return -q[index3_to_index(uint3(index.x,index.y,0))];
else if(index.z >= FLUID_GRID_CELL_COUNT.z)
return -q[index3_to_index(uint3(index.x,index.y,FLUID_GRID_CELL_COUNT.z-1))];
return q[index3_to_index(index)];
}
[numthreads(1, 1, 1)]
void main( uint3 DTid : SV_DispatchThreadID )
{
int cell_index_center = index3_to_index(DTid);
float3 q_center = q[cell_index_center];
float3 q_left = get_quantity_bd(DTid + uint3(-1,0,0));
float3 q_right = get_quantity_bd(DTid + uint3(1,0,0));
float3 q_bottom = get_quantity_bd(DTid + uint3(0,-1,0));
float3 q_top = get_quantity_bd(DTid + uint3(0,1,0));
float3 q_front = get_quantity_bd(DTid + uint3(0,0,-1));
float3 q_back = get_quantity_bd(DTid + uint3(0,0,1));
float div = 0.5 * FLUID_GRID_CELL_CONST.x * ((q_right.x - q_left.x)+
(q_top.y - q_bottom.y)+
(q_back.z - q_front.z));
divergence[cell_index_center] = div;
}
Code: Alles auswählen
float get_quantity_bd(uint3 index)
{
float value;
float3 mod_index;
if(index.x < 0)
mod_index.x = 0;
else if(index.x >= FLUID_GRID_CELL_COUNT.x)
mod_index.x = FLUID_GRID_CELL_COUNT.x-1;
else
mod_index.x = index.x;
if(index.y < 0)
mod_index.y = 0;
else if(index.y >= FLUID_GRID_CELL_COUNT.y)
mod_index.y = FLUID_GRID_CELL_COUNT.y-1;
else
mod_index.y = index.y;
if(index.z < 0)
mod_index.z = 0;
else if(index.z >= FLUID_GRID_CELL_COUNT.z)
mod_index.z = FLUID_GRID_CELL_COUNT.z-1;
else
mod_index.z = index.z;
value = q[index3_to_index(mod_index)];
return value;
}
[numthreads(1, 1, 1)]
void main( uint3 DTid : SV_DispatchThreadID )
{
int cell_index_center = index3_to_index(DTid); //OK
float div = b[cell_index_center];
float q_center = q[cell_index_center];
float q_left = get_quantity_bd(DTid + uint3(-1,0,0));
float q_right = get_quantity_bd(DTid + uint3(1,0,0));
float q_bottom = get_quantity_bd(DTid + uint3(0,-1,0));
float q_top = get_quantity_bd(DTid + uint3(0,1,0));
float q_front = get_quantity_bd(DTid + uint3(0,0,1));
float q_back = get_quantity_bd(DTid + uint3(0,0,-1));
x[cell_index_center] = (q_left + q_right + q_bottom + q_top + q_front + q_back - div) / 6.0;
}
Code: Alles auswählen
float get_quantity_bd(uint3 index)
{
float value;
float3 mod_index;
if(index.x < 0)
mod_index.x = 0;
else if(index.x >= FLUID_GRID_CELL_COUNT.x)
mod_index.x = FLUID_GRID_CELL_COUNT.x-1;
else
mod_index.x = index.x;
if(index.y < 0)
mod_index.y = 0;
else if(index.y >= FLUID_GRID_CELL_COUNT.y)
mod_index.y = FLUID_GRID_CELL_COUNT.y-1;
else
mod_index.y = index.y;
if(index.z < 0)
mod_index.z = 0;
else if(index.z >= FLUID_GRID_CELL_COUNT.z)
mod_index.z = FLUID_GRID_CELL_COUNT.z-1;
else
mod_index.z = index.z;
value = q[index3_to_index(mod_index)];
return value;
}
[numthreads(1, 1, 1)]
void main( uint3 DTid : SV_DispatchThreadID )
{
int cell_index_center = index3_to_index(DTid);
float q_left = get_quantity_bd(DTid + uint3(-1,0,0));
float q_right = get_quantity_bd(DTid + uint3(1,0,0));
float q_bottom = get_quantity_bd(DTid + uint3(0,-1,0));
float q_top = get_quantity_bd(DTid + uint3(0,1,0));
float q_front = get_quantity_bd(DTid + uint3(0,0,-1));
float q_back = get_quantity_bd(DTid + uint3(0,0,1));
float3 gradP = 0.5*float3(q_right - q_left, q_top - q_bottom, q_back - q_front);
velocity[cell_index_center] = velocity_old[cell_index_center] - gradP;
}
Wenn ich nach der Projektion das Div Feld neu berechnen ist leider nicht wirklich divergenzfrei.
Vielen Dank,
Manuel