Ocean Simulation

Unity / C# & HLSL

Physics & Simulation

What Is Happening Here

There are two main parts to this project. One is the ocean that is comprised of waves moving along the surface and the other is the physics object that inhabits this ocean, in this case a small beach ball. The goal of the simulation is to recreate physics like behaviour in a dynamic ocean we have some parts down and here I will describe what I did to solve the problems I have encountered on the way.


Waves

So what is a wave, Webster’s dictionary says: “to motion with the hands or with something held in them in signal or salute”. Well that is not quite what I am looking for, what we want is an oscillating movement, the same thing we have in things like soundwaves and pendulums. So, what happens if we create a basic sine curve, we get something like the gif on the right.

But normally waves don’t look as uniform as this, they have more of a crest and steepness to them. So, we somehow have to make them act more naturally. This is where the Gerstner-wave comes in, this is what is represented on the graph. When adding two sine waves you can create more complex characteristics of the waves. In my example I have only implemented a basic form of Gerster-waves to give it some steepness in the future I am looking to add more directions to give it a more realistic feeling with different sizes of waves and patterns, at this point we have this uniform wave pattern. Now that we know what sort of pattern we want our wave to create, how do we go about and create this in code? First of all we should start of with generating of a plane so that we have some vertices to manipulate in to our decried pattern. Remember we want this to be able to execute in runtime for us to be able to use this data for physics calculations. It just so happens to be that we have a gpu that is quite good at handling vast amount of float values. Here is how I did it

[numthreads(1,1,1)]
void GenerateMesh(uint3 id : SV_DispatchThreadID)
{
    float offset = float(SeaResolution * SeaSize * 0.5f);
    float maxSize = SeaSize * SeaResolution;
	//Creating verticies & Uvs
	for (uint i = 0; i < VerticesPerSide; i++)
	{
		for (uint j = 0; j < VerticesPerSide; j++)
		{
			Vertices[i * VerticesPerSide + j] = float3(j * SeaResolution, 0, i * SeaResolution) - offset;
			Uvs[i * VerticesPerSide + j] = (float2(j * SeaResolution, i * SeaResolution) - offset);
		}
	}

	int nrVerticies = VerticesPerSide * VerticesPerSide;

	uint2 index = uint2(0,0);
	uint triangleIndex = 0;
	int size = nrVerticies * 6;

	//Creating triangles
	for (int k = 0; k < nrVerticies; k++)
	{

		index.x = k / VerticesPerSide;
		index.y = k % VerticesPerSide;

		if (index.x == SeaSize || index.y == SeaSize)
			continue;

		Triangles[triangleIndex++] = index.x + index.y * VerticesPerSide; // (0, 0) bottom left
		Triangles[triangleIndex++] = index.x + (index.y + 1) * VerticesPerSide; // (0, 1) top left
		Triangles[triangleIndex++] = index.x + 1 + (index.y + 1) * VerticesPerSide; // (1, 1) top right

		Triangles[triangleIndex++] = index.x + 1 + (index.y + 1) * VerticesPerSide; // (1, 1) top right
		Triangles[triangleIndex++] = index.x + 1 + index.y * VerticesPerSide; // (1, 0) bottom right
		Triangles[triangleIndex++] = index.x + index.y * VerticesPerSide; // (0, 0) bottom left
	}
}
                

So now we have our plane lets move those vertices! But moving them on the cpu will be very slow and we can’t just do it in a regular shader because we want to access the vertices height and position for our calculations. Well there just happens to be a solution, compute shaders. A compute shader in and of itself does not handle meshes but if we feed it the mesh data it can do the calculations and give it back to the cpu. The cpu tells the gpu what the current time of the simulation is for it to temporally change and what the vertices it is working on. Then it is only a matter of crunching the numbers.


So now we have our plane lets move those vertices! But moving them on the cpu will be very slow and we can’t just do it in a regular shader because we want to access the vertices height and position for our calculations. Well there just happens to be a solution, compute shaders. A compute shader in and of itself does not handle meshes but if we feed it the mesh data it can do the calculations and give it back to the cpu. The cpu tells the gpu what the current time of the simulation is for it to temporally change and what the vertices it is working on. Then it is only a matter of crunching the numbers.

RWStructuredBuffer ActiveVertices;
float WaveLength = 0.5;
float Amplitude;
float Speed;
float2 Direction;
uint SizePerSide;
float Resolution;
float Frequency;

float Time;
float Steepness;
float Gerst;

[numthreads(1, 1, 1)]
void UpdateMesh(uint3 id : SV_DispatchThreadID)
{
	float3 currentVector  = ActiveVertices[id.x + id.y * SizePerSide];
	float2 undisturbedVector  = float2(id.x,  id.y);
   // float b = 1 - (Steepness * dot(normalize(Direction)
    float k = 6.28318530718f / WaveLength;
    float w = dot(normalize(Direction) * k, undisturbedVector) - (6.28318530718f / Frequency) * Time;
    currentVector.x = undisturbedVector.x - (normalize(Direction) * Amplitude * sin(w));
    currentVector.y = Amplitude * cos(w);
    //
    currentVector.x -= Resolution * SizePerSide * 0.5f;

    ActiveVertices[id.x + id.y * SizePerSide] = currentVector;
}
                

Buoyancy

Now that we have some ocean looking pattern moving along our plane it is time for us to put those values that we received from the gpu to use. The formula for calculating the buoyance is quite trivial if we have a sphere, that is why we will be simulating a beach ball. All we need is the radius and density of the ball and how submerged it is to create buoyancy force. This is where the height values form the vertices come in. We can calculate where on the ocean the ball is and get the current height of the ocean at the point. Here is how I simulated buoyancy in my project.

public class Buoyancy : MonoBehaviour
{
    [SerializeField]private float _radius = 1;
    [SerializeField]private float _density = 1;
    private float _depth;
    private static readonly float RHO_GRAVITY = Physics.gravity.magnitude * 1000f;  //rho = density of water = 1000 Kg/m3 , gravity = 10 m/s2
    private Rigidbody _rigidbody;
    private const float _densityOfAir = 1.2754f;
    private const float _densityOfWater = 997f;
    private const float _sphereDragCoefficient = 0.47f;
    private Plane _plane;
    void Start()
    {
        _rigidbody = gameObject.AddComponent();
        _rigidbody.mass = GetVolume(_radius) * _density;
        _plane = new Plane();
    }

    void FixedUpdate()
    {
        float buoyancyForce = RHO_GRAVITY * GetSubmersedVolume(_radius);
        CalculateDrag();
        //Debug.Log(buoyancyForce);
        _rigidbody.AddForce(Vector3.up * buoyancyForce, ForceMode.Force);

    }

    private float GetSubmersedVolume(float radius)
    {
        float p = Vector3.Dot(SeaController.GetSurfaceNormal(transform.position),transform.position - SeaController.GetVertexPosition(transform.position));
        _depth = radius - p;
        _depth = Mathf.Clamp(_depth, 0f, radius * 2);
        return Mathf.PI * ((radius * (_depth * _depth)) - (_depth * _depth * _depth) / 3f) ;
    }

    private void CalculateDrag()
    {

        //Do air drag
        if ((transform.position.y - _radius )- SeaController.GetHeight(transform.position) > 0)
        {
            float dragForce = (_densityOfAir / 2) * _rigidbody.velocity.magnitude * _rigidbody.velocity.magnitude * _sphereDragCoefficient * _radius * _radius * Mathf.PI;
            _rigidbody.AddForce(dragForce * (-_rigidbody.velocity.normalized));
        }
        //Do Water drag
        else
        {
            float dragForce = (_densityOfWater / 2) * _rigidbody.velocity.magnitude * _rigidbody.velocity.magnitude * _sphereDragCoefficient * _radius * _radius * Mathf.PI;
            _rigidbody.AddForce(dragForce * (-_rigidbody.velocity.normalized));
        }
    }

    private float GetVolume(float radius)
    {
        return (4f/3f) * Mathf.PI * radius * radius * radius;
    }

    private float GetVolume(Mesh mesh)
    {
        float volume = 0;

        float lowestY = 0;


        return volume;
    }
}