Can you fix this fluid simulation.  Currently it just makes one big column of water and then goes crazy. I want it to flow on the heightmap of the terrain. Here is the code: --[[
    Height Field Fluid Simulation in Roblox Studio
    Implementing the methods from "Height Field Fluids" by Matthias Müller
]]

-- Grid dimensions
local GRID_WIDTH = 100    -- Adjust as needed
local GRID_HEIGHT = 100   -- Adjust as needed

-- Grid spacing (size of each grid cell in studs)
local DX = 4              -- Adjust based on desired resolution
local DY = 4

-- Physical constants
local GRAVITY = 9.81      -- Acceleration due to gravity (studs/s^2)

-- Simulation parameters
local TIME_STEP = 0.03    -- Initial time step for the simulation (in seconds)
local MIN_TIME_STEP = 0.001 -- Minimum allowable time step to prevent zero division
local MAX_TIME_STEP = 0.03  -- Maximum allowable time step
local CFL_NUMBER = 0.5    -- Courant number for CFL condition (usually <= 1)

-- Height fields and velocities
local H = {}           -- Total height (fluid + terrain)
local H_Terrain = {}   -- Terrain height
local H_Water = {}     -- Water depth (fluid height)
local U = {}           -- Momentum in x-direction
local V = {}           -- Momentum in y-direction

-- Fluxes
local Flux_X = {}      -- Flux in x-direction
local Flux_Y = {}      -- Flux in y-direction

-- Initialize the data structures
for i = 1, GRID_WIDTH + 1 do
	H[i] = {}
	H_Terrain[i] = {}
	H_Water[i] = {}
	U[i] = {}
	V[i] = {}
	Flux_X[i] = {}
	Flux_Y[i] = {}
	for j = 1, GRID_HEIGHT + 1 do
		H[i][j] = 0
		H_Terrain[i][j] = 0
		H_Water[i][j] = 0
		U[i][j] = 0
		V[i][j] = 0
		Flux_X[i][j] = { h = 0, hu = 0, hv = 0 }
		Flux_Y[i][j] = { h = 0, hu = 0, hv = 0 }
	end
end

-- Function to sample terrain height at a given position using raycasts
local function SampleTerrainHeight(x, y)
	-- Convert grid coordinates to world coordinates
	local worldX = (x - 0.5) * DX
	local worldZ = (y - 0.5) * DY
	local origin = Vector3.new(worldX, 500, worldZ)  -- Start raycast from above (500 studs high)
	local direction = Vector3.new(0, -1000, 0)       -- Raycast downwards (1000 studs)

	-- Create a raycast parameters object
	local params = RaycastParams.new()
	params.FilterDescendantsInstances = {}  -- Add instances to ignore if needed
	params.FilterType = Enum.RaycastFilterType.Blacklist
	params.IgnoreWater = false  -- Include water in raycast if needed

	local rayResult = workspace:Raycast(origin, direction, params)

	if rayResult then
		local hitPosition = rayResult.Position
		local terrainHeight = hitPosition.Y
		return terrainHeight
	else
		-- Default terrain height if nothing is hit
		return 0
	end
end

-- Generate the terrain heightmap
function GenerateTerrainHeightmap()
	for i = 1, GRID_WIDTH + 1 do
		for j = 1, GRID_HEIGHT + 1 do
			H_Terrain[i][j] = SampleTerrainHeight(i, j)
			H_Water[i][j] = 0
			H[i][j] = H_Terrain[i][j]
			U[i][j] = 0
			V[i][j] = 0
		end
	end
end

-- Water Source
local WaterSource = {
	x = math.floor(GRID_WIDTH / 2),
	y = math.floor(GRID_HEIGHT / 2),
	rate = 1  -- Water addition rate (studs per second)
}

-- Add water at the source
function AddWater(deltaTime)
	local i, j = WaterSource.x, WaterSource.y
	local rate = WaterSource.rate * deltaTime
	H_Water[i][j] = H_Water[i][j] + rate
end

-- Compute maximum wave speed for CFL condition
function ComputeMaxWaveSpeed()
	local maxSpeed = 0
	for i = 2, GRID_WIDTH do
		for j = 2, GRID_HEIGHT do
			local h = H_Water[i][j]
			if h > 0 then
				local u = U[i][j] / h
				local v = V[i][j] / h
				local speed = math.sqrt(u * u + v * v) + math.sqrt(GRAVITY * h)
				if speed > maxSpeed then
					maxSpeed = speed
				end
			end
		end
	end
	return maxSpeed
end

-- Compute time step based on CFL condition
function ComputeTimeStep()
	local maxSpeed = ComputeMaxWaveSpeed()
	if maxSpeed > 0 then
		local dt = CFL_NUMBER * (math.min(DX, DY) / maxSpeed)
		dt = math.clamp(dt, MIN_TIME_STEP, MAX_TIME_STEP)
		return dt
	else
		return TIME_STEP  -- Default time step if no movement
	end
end

-- Compute fluxes using HLL approximate Riemann solver
function ComputeFluxes()
	-- Compute fluxes in x-direction
	for i = 1, GRID_WIDTH do
		for j = 1, GRID_HEIGHT do
			-- Left (L) and right (R) states
			local hL = H_Water[i][j]
			local hR = H_Water[i + 1][j]
			local uL = (hL > 0) and (U[i][j] / hL) or 0
			local uR = (hR > 0) and (U[i + 1][j] / hR) or 0
			local vL = (hL > 0) and (V[i][j] / hL) or 0
			local vR = (hR > 0) and (V[i + 1][j] / hR) or 0

			-- Compute wave speeds
			local sqrt_g_hL = math.sqrt(GRAVITY * hL)
			local sqrt_g_hR = math.sqrt(GRAVITY * hR)
			local sL = math.min(uL - sqrt_g_hL, uR - sqrt_g_hR)
			local sR = math.max(uL + sqrt_g_hL, uR + sqrt_g_hR)

			-- Fluxes for h, hu, hv
			local fluxL = {
				h = U[i][j],
				hu = U[i][j] * uL + 0.5 * GRAVITY * hL * hL,
				hv = U[i][j] * vL
			}
			local fluxR = {
				h = U[i + 1][j],
				hu = U[i + 1][j] * uR + 0.5 * GRAVITY * hR * hR,
				hv = U[i + 1][j] * vR
			}

			-- HLL Flux
			if sL >= 0 then
				Flux_X[i][j] = fluxL
			elseif sR <= 0 then
				Flux_X[i][j] = fluxR
			else
				local denom = sR - sL
				Flux_X[i][j] = {
					h = (sR * fluxL.h - sL * fluxR.h + sL * sR * (hR - hL)) / denom,
					hu = (sR * fluxL.hu - sL * fluxR.hu + sL * sR * (U[i + 1][j] - U[i][j])) / denom,
					hv = (sR * fluxL.hv - sL * fluxR.hv + sL * sR * (V[i + 1][j] - V[i][j])) / denom
				}
			end
		end
	end

	-- Compute fluxes in y-direction
	for i = 1, GRID_WIDTH do
		for j = 1, GRID_HEIGHT do
			-- Bottom (B) and top (T) states
			local hB = H_Water[i][j]
			local hT = H_Water[i][j + 1]
			local uB = (hB > 0) and (U[i][j] / hB) or 0
			local uT = (hT > 0) and (U[i][j + 1] / hT) or 0
			local vB = (hB > 0) and (V[i][j] / hB) or 0
			local vT = (hT > 0) and (V[i][j + 1] / hT) or 0

			-- Compute wave speeds
			local sqrt_g_hB = math.sqrt(GRAVITY * hB)
			local sqrt_g_hT = math.sqrt(GRAVITY * hT)
			local sB = math.min(vB - sqrt_g_hB, vT - sqrt_g_hT)
			local sT = math.max(vB + sqrt_g_hB, vT + sqrt_g_hT)

			-- Fluxes for h, hu, hv
			local fluxB = {
				h = V[i][j],
				hu = V[i][j] * uB,
				hv = V[i][j] * vB + 0.5 * GRAVITY * hB * hB
			}
			local fluxT = {
				h = V[i][j + 1],
				hu = V[i][j + 1] * uT,
				hv = V[i][j + 1] * vT + 0.5 * GRAVITY * hT * hT
			}

			-- HLL Flux
			if sB >= 0 then
				Flux_Y[i][j] = fluxB
			elseif sT <= 0 then
				Flux_Y[i][j] = fluxT
			else
				local denom = sT - sB
				Flux_Y[i][j] = {
					h = (sT * fluxB.h - sB * fluxT.h + sB * sT * (hT - hB)) / denom,
					hu = (sT * fluxB.hu - sB * fluxT.hu + sB * sT * (U[i][j + 1] - U[i][j])) / denom,
					hv = (sT * fluxB.hv - sB * fluxT.hv + sB * sT * (V[i][j + 1] - V[i][j])) / denom
				}
			end
		end
	end
end

-- Update conserved variables (water depth and momenta)
function UpdateConservedVariables(deltaTime)
	for i = 2, GRID_WIDTH do
		for j = 2, GRID_HEIGHT do
			-- Net fluxes
			local flux_h = (Flux_X[i - 1][j].h - Flux_X[i][j].h) / DX + (Flux_Y[i][j - 1].h - Flux_Y[i][j].h) / DY
			local flux_hu = (Flux_X[i - 1][j].hu - Flux_X[i][j].hu) / DX + (Flux_Y[i][j - 1].hu - Flux_Y[i][j].hu) / DY
			local flux_hv = (Flux_X[i - 1][j].hv - Flux_X[i][j].hv) / DX + (Flux_Y[i][j - 1].hv - Flux_Y[i][j].hv) / DY

			-- Update water depth
			H_Water[i][j] = H_Water[i][j] - deltaTime * flux_h

			-- Update momenta
			U[i][j] = U[i][j] - deltaTime * flux_hu
			V[i][j] = V[i][j] - deltaTime * flux_hv

			-- Ensure non-negative water depth
			if H_Water[i][j] < 0 then
				H_Water[i][j] = 0
				U[i][j] = 0
				V[i][j] = 0
			end
		end
	end
end

-- Enforce boundary conditions
function EnforceBoundaryConditions()
	-- Reflective boundary conditions for left and right boundaries
	for j = 1, GRID_HEIGHT + 1 do
		H_Water[1][j] = H_Water[2][j]
		U[1][j] = -U[2][j]
		V[1][j] = V[2][j]

		H_Water[GRID_WIDTH + 1][j] = H_Water[GRID_WIDTH][j]
		U[GRID_WIDTH + 1][j] = -U[GRID_WIDTH][j]
		V[GRID_WIDTH + 1][j] = V[GRID_WIDTH][j]
	end

	-- Reflective boundary conditions for top and bottom boundaries
	for i = 1, GRID_WIDTH + 1 do
		H_Water[i][1] = H_Water[i][2]
		U[i][1] = U[i][2]
		V[i][1] = -V[i][2]

		H_Water[i][GRID_HEIGHT + 1] = H_Water[i][GRID_HEIGHT]
		U[i][GRID_HEIGHT + 1] = U[i][GRID_HEIGHT]
		V[i][GRID_HEIGHT + 1] = -V[i][GRID_HEIGHT]
	end
end

-- Update total height (for visualization)
function UpdateTotalHeight()
	for i = 1, GRID_WIDTH + 1 do
		for j = 1, GRID_HEIGHT + 1 do
			H[i][j] = H_Terrain[i][j] + H_Water[i][j]
		end
	end
end

-- Rendering the simulation
local WaterParts = {}

function InitializeWaterParts()
	local WaterModel = Instance.new("Model")
	WaterModel.Name = "Water"
	WaterModel.Parent = workspace

	for i = 1, GRID_WIDTH do
		WaterParts[i] = {}
		for j = 1, GRID_HEIGHT do
			local part = Instance.new("Part")
			part.Size = Vector3.new(DX, 1, DY)
			part.Anchored = true
			part.CanCollide = false
			part.Material = Enum.Material.Water
			part.Color = Color3.fromRGB(0, 162, 255)
			part.Transparency = 0.5
			part.Parent = WaterModel
			WaterParts[i][j] = part
		end
	end
end

function RenderSimulation()
	for i = 1, GRID_WIDTH do
		for j = 1, GRID_HEIGHT do
			local waterDepth = H_Water[i][j]
			local terrainHeight = H_Terrain[i][j]
			local part = WaterParts[i][j]

			if waterDepth > 0.1 then
				-- Position the water part
				part.Position = Vector3.new(
					(i - 0.5) * DX,
					terrainHeight + waterDepth / 2,
					(j - 0.5) * DY
				)
				-- Adjust the size (height)
				part.Size = Vector3.new(DX, waterDepth, DY)
				part.Transparency = 0.5

				-- Optionally, visualize flow velocity using color
				local h = H_Water[i][j]
				local u = (h > 0) and (U[i][j] / h) or 0
				local v = (h > 0) and (V[i][j] / h) or 0
				local speed = math.sqrt(u * u + v * v)
				local colorIntensity = math.clamp(speed / 20, 0, 1)  -- Adjust divisor as needed
				part.Color = Color3.new(0, 0.5 + 0.5 * colorIntensity, 1)
			else
				-- Hide the part if no water
				part.Size = Vector3.new(0, 0, 0)
				part.Transparency = 1
			end
		end
	end
end

-- Initialize the simulation
GenerateTerrainHeightmap()
InitializeWaterParts()

-- Main simulation loop
while true do
	local deltaTime = ComputeTimeStep()

	-- Add water at the source
	AddWater(deltaTime)

	-- Compute fluxes using HLL Riemann solver
	ComputeFluxes()

	-- Update conserved variables (H_Water, U, V)
	UpdateConservedVariables(deltaTime)

	-- Enforce boundary conditions
	EnforceBoundaryConditions()

	-- Update total height (for visualization)
	UpdateTotalHeight()

	-- Render the simulation
	RenderSimulation()

	-- Wait for the next time step
	wait(deltaTime)
end