Geschwindichkeitsfeld divergenz frei machen

Design Patterns, Erklärungen zu Algorithmen, Optimierung, Softwarearchitektur
Forumsregeln
Wenn das Problem mit einer Programmiersprache direkt zusammenhängt, bitte HIER posten.
Antworten
MadMax
Beiträge: 59
Registriert: 24.01.2003, 13:31
Kontaktdaten:

Geschwindichkeitsfeld divergenz frei machen

Beitrag von MadMax »

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:

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;
	
}
2.) Jacob Iteration (x500):

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;  
	
}
3.) Projektion:

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;  		
	
}
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
Zuletzt geändert von MadMax am 05.08.2012, 12:20, insgesamt 2-mal geändert.
MadMax
Beiträge: 59
Registriert: 24.01.2003, 13:31
Kontaktdaten:

Re: Geschwindichkeitsfeld divergenz frei machen

Beitrag von MadMax »

Fehler Nr.1 gerade gefunden:
Alle uint durch int ersetzt....

Funktioniert aber trotdem nicht....
MadMax
Beiträge: 59
Registriert: 24.01.2003, 13:31
Kontaktdaten:

Re: Geschwindichkeitsfeld divergenz frei machen

Beitrag von MadMax »

Fehler Nr. 2:

die funktion die die Randwerte berechnet ist falsch gewesen.

Code: Alles auswählen

float3 get_quantity_bd(int3 index)
{
	float3 value;
	
	if(index.x < 0)
		return float3(-q[index3_to_index(int3(0,index.y,index.z))].x,0.0,0.0); 
	
	if(index.x >= FLUID_GRID_CELL_COUNT.x)
		return float3(-q[index3_to_index(int3(FLUID_GRID_CELL_COUNT.x-1,index.y,index.z))].x,0.0,0.0);

	if(index.y < 0)
		return float3(0.0,-q[index3_to_index(int3(index.x,0,index.z))].y,0.0); 
	
	if(index.y >= FLUID_GRID_CELL_COUNT.y)
		return float3(0.0,-q[index3_to_index(int3(index.x,FLUID_GRID_CELL_COUNT.y-1,index.z))].y,0.0);

	if(index.z < 0)
		return float3(0.0,0.0,-q[index3_to_index(int3(index.x,index.y,0))].z); 
	
	if(index.z >= FLUID_GRID_CELL_COUNT.y)
		return float3(0.0,0.0,-q[index3_to_index(int3(index.x,index.y,FLUID_GRID_CELL_COUNT.z-1))].z);

	return q[index3_to_index(index)];
}
Funktioniert leider immer noch nicht...
Zuletzt geändert von MadMax am 05.08.2012, 12:39, insgesamt 2-mal geändert.
MadMax
Beiträge: 59
Registriert: 24.01.2003, 13:31
Kontaktdaten:

Re: Geschwindichkeitsfeld divergenz frei machen

Beitrag von MadMax »

Schade hat sich noch niemand mit Flüßigkeitssimulationen beschäftigt. Ich bin über jede Anregung dankbar.
Benutzeravatar
Artificial Mind
Establishment
Beiträge: 802
Registriert: 17.12.2007, 17:51
Wohnort: Aachen

Re: Geschwindichkeitsfeld divergenz frei machen

Beitrag von Artificial Mind »

Naja, beschäftigt schon. Meine Bachelorarbeit schreibe ich gerade über explizite Integration der Navier-Stokes-Gleichungen für den kompressiblen Fall. Da ich hier aber nicht mehr die Inkompressibilitätsannahme habe, habe ich auch kein divergenzfreies Geschindigkeitsfeld.

Ich finde deinen Code äußerst anstrengend zu lesen (kein Highlighting, keine Einrückungen, minimalistische Variablennamen, keine Kommentare), deswegen konnte ich mich noch nicht überwinden, mir das genau durchzulesen.

Das Einzige, was ich dir sagen kann, ist, dass du wahrscheinlich die Stable Fluids versuchst, nachzuprogrammieren. Dort kenne ich diese Webseite: http://www.multires.caltech.edu/teachin ... fluids.htm
Ist zwar in Java, aber ich hatte den Code auch schon mal nach C++ portiert und das funktioniert gut.
MadMax
Beiträge: 59
Registriert: 24.01.2003, 13:31
Kontaktdaten:

Re: Geschwindichkeitsfeld divergenz frei machen

Beitrag von MadMax »

Was Einrückungen und Highlighting angeht gebe ich dir recht. Aber das haben bei mir die code tags gemacht. Ka wie ich das verhindern kann.
Benutzeravatar
Artificial Mind
Establishment
Beiträge: 802
Registriert: 17.12.2007, 17:51
Wohnort: Aachen

Re: Geschwindichkeitsfeld divergenz frei machen

Beitrag von Artificial Mind »

statt nur "code" kannst du "code=cpp" oder "code=hlsl" oder so versuchen. Irgendwelche davon produzieren wesentlich besseres Highlighting.
Antworten