1- package io .computenode .cyfra .fluids .solver
1+ package io .computenode .cyfra .fluids .solver . programs
22
33import io .computenode .cyfra .core .GProgram
44import io .computenode .cyfra .core .GProgram .StaticDispatch
55import io .computenode .cyfra .dsl .{* , given }
6- import io .computenode .cyfra .fluids .solver .GridUtils .idxTo3D
6+ import io .computenode .cyfra .fluids .solver .*
7+ import io .computenode .cyfra .fluids .solver .utils .{GridUtils , ObstacleUtils }
8+ import GridUtils .idxTo3D
79
810/** Free-slip boundary conditions with density/temperature dissipation.
911 *
@@ -46,7 +48,6 @@ object DissipativeBoundaryProgram:
4648 (z === 0 ) || (z === n - 1 )
4749 val isSolid = ObstacleUtils .isSolid(state.obstacles, idx, totalCells)
4850
49- // Check if adjacent to an obstacle (for obstacle surface boundary conditions)
5051 val solidXP = ObstacleUtils .isSolidAt(state.obstacles, x + 1 , y, z, n)
5152 val solidXM = ObstacleUtils .isSolidAt(state.obstacles, x - 1 , y, z, n)
5253 val solidYP = ObstacleUtils .isSolidAt(state.obstacles, x, y + 1 , z, n)
@@ -55,59 +56,48 @@ object DissipativeBoundaryProgram:
5556 val solidZM = ObstacleUtils .isSolidAt(state.obstacles, x, y, z - 1 , n)
5657 val adjacentToObstacle = (solidXP || solidXM || solidYP || solidYM || solidZP || solidZM) && ! isSolid
5758
58- // Compose both GIO operations
5959 for
60- // 1. Handle obstacles - zero out everything inside
6160 _ <- GIO .when(isSolid):
6261 for
6362 _ <- GIO .write(state.velocity, idx, vec4(0.0f , 0.0f , 0.0f , 0.0f ))
6463 _ <- GIO .write(state.density, idx, 0.0f )
6564 _ <- GIO .write(state.temperature, idx, 0.0f )
6665 yield GStruct .Empty ()
6766
68- // 2. Free-slip at obstacle surfaces (fluid cells adjacent to obstacles)
6967 _ <- GIO .when(adjacentToObstacle && ! onDomainBoundary):
7068 val vel = state.velocity.read(idx)
7169 val vel3 = vec3(vel.x, vel.y, vel.z)
7270
73- // Compute normal pointing away from obstacle (sum of normals for all adjacent solid faces)
7471 val nx = when(solidXP)(1.0f ).otherwise(0.0f ) + when(solidXM)(- 1.0f ).otherwise(0.0f )
7572 val ny = when(solidYP)(1.0f ).otherwise(0.0f ) + when(solidYM)(- 1.0f ).otherwise(0.0f )
7673 val nz = when(solidZP)(1.0f ).otherwise(0.0f ) + when(solidZM)(- 1.0f ).otherwise(0.0f )
7774 val obstacleNormal = vec3(nx, ny, nz)
7875
79- // Normalize
8076 val normalLength = sqrt((obstacleNormal dot obstacleNormal) + 1e-8f )
8177 val normalNorm = obstacleNormal * (1.0f / normalLength)
8278
83- // Free-slip: remove normal component
8479 val normalComponent = vel3 dot normalNorm
8580 val tangentialVel = vel3 - (normalNorm * normalComponent)
8681 val newVel = vec4(tangentialVel.x, tangentialVel.y, tangentialVel.z, 0.0f )
8782
8883 GIO .write(state.velocity, idx, newVel)
8984
90- // 3. Free-slip with dissipation at domain boundaries
9185 _ <- GIO .when(onDomainBoundary && ! isSolid && ! adjacentToObstacle):
9286 val vel = state.velocity.read(idx)
9387 val vel3 = vec3(vel.x, vel.y, vel.z)
9488
95- // Domain boundary normals
9689 val nx = when(x === 0 )(1.0f ).elseWhen(x === n - 1 )(- 1.0f ).otherwise(0.0f )
9790 val ny = when(y === 0 )(1.0f ).elseWhen(y === n - 1 )(- 1.0f ).otherwise(0.0f )
9891 val nz = when(z === 0 )(1.0f ).elseWhen(z === n - 1 )(- 1.0f ).otherwise(0.0f )
9992 val normal = vec3(nx, ny, nz)
10093
101- // Normalize for corners (where multiple components are non-zero)
10294 val normalLength = sqrt((normal dot normal) + 1e-8f )
10395 val normalNorm = normal * (1.0f / normalLength)
10496
105- // Free-slip: remove normal component
10697 val normalComponent = vel3 dot normalNorm
10798 val tangentialVel = vel3 - (normalNorm * normalComponent)
10899 val newVel = vec4(tangentialVel.x, tangentialVel.y, tangentialVel.z, 0.0f )
109100
110- // Dissipate density and temperature at domain boundaries
111101 val currentDensity = state.density.read(idx)
112102 val currentTemp = state.temperature.read(idx)
113103 val dissipationFactor = 0.95f
0 commit comments