January 8, 2018 at 21:26GPU physical simulation using compute shader in Unity3D
- Tutorial
In this tutorial, I’ll show you how to use compute shader to implement computations on a video card - using an example of a hair model:
Here is a project for Unity3D, on the explanation of which the guide is built. You need to download and open it in Unity: a
link to the Unity project.
Who will understand this guide? Those who use Unity3D or at least know C # or C ++. The shader is written in HLSL, a close syntactic cousin of C ++.
Who will benefit from this guide? Experienced programmers who want to learn how to use the GPU for computing. But even an inexperienced but diligent programmer will easily understand everything.
Why use a graphics card for computing? For parallel tasks, its performance is 10-100 times higher than that of the processor. That is, everyone in the computer has a small supercomputer with a convenient API, it makes sense to use it in suitable cases.
Is this huge performance really needed? Yes, processor speed is often a limiting factor. For example, when you need to perform the same operations on large data sets. But it is precisely such tasks that can be easily parallelized. In addition, developers often refuse decisions because of their computing capacity, and entire areas in the algorithm space remain unexplored. For example, you can do the coolest physics in games, if you properly load the GPU.
And what, with a video card you can now simply solve problems with brute force? The demand for optimization does not depend on the performance of iron. There is no such supercomputer that could not be tightly loaded with inefficient code.
Why exactly compute shader? Why not opencl or cuda? Cuda only works on nvidia hardware, and I don’t know opencl. Unity can build in any API, including opengl core. On Macs and on an Android computer, shaders work, on Linux, it seems, too (although I have not tried it). Although, each API has limitations that should be considered. For example, on Metal you cannot make more than 256 threads along one axis (In DX10 - 1024). And the android API will not be able to use more than 4 buffers per kernel (in DX10 - 8, in DX11 - even more).
Why physical simulation? This is a computationally intensive task, and is well suited for parallel computing. In addition, the task is in demand. Game devs can implement interesting physics in games, students can create experimental models for term papers, engineers and scientists can do calculations on models.
And why exactly the hair model? I wanted to take a simple puzzle, but at the same time covering the main problems.
How to use this manual? It is best to download the source code, open it and read it as you progress through the manual. I will explain in detail all the main lines, although I will not explain absolutely every line, the meaning of most of them is obvious. There are no complex algorithms in the text, there is only the use of the interface of classes that serve computing on the GPU. And on the side of the shader code, there is nothing but reading the data, performing simple mathematical operations on them, and recording the results. But if something is not clear - be sure to ask, I will answer everything in kamenty.
And now for those who have absolutely no idea about using compute shaders, I suggest taking a step to the side and moving on to a very simple guide, which is devoted to the basics of using computer shaders. I advise you to start with it in order to better understand the essence and get used to the practice of GPU computing with an extremely simple example. And then come back here and continue. And those who are at least somehow familiar with computer shaders, let them boldly read on.
If you want to make a physical model calculated on a GPU from scratch, then this problem can be divided into 4 parts:
- a mathematical model of the phenomenon
- an algorithm for parallel model calculation
- a shader code
- preparing and running a shader in a unit
The strengths of video cards are that they can apply one operation simultaneously to many objects. Therefore, a hair model can be made as a chain of dots, each of which interacts with two neighbors. The interaction between the points is based on the spring principle: k * (S0-S) ^ n, where S0 is the equilibrium distance, S is the current distance. In reality, the hair does not look like a spring, it is perceived as inextensible. This means that the spring in the model must be made sufficiently rigid. It is better to increase the stiffness of the spring by increasing n, because the degree increases the curvature of the curve in the vicinity of the equilibrium, which reduces the backlash and reduces the effect of “rubberiness” of the hair. I took n = 2, and we will talk about the value of the coefficient k below.
In addition to the elastic force between the points, diffusion of relative velocities or one-dimensional viscosity will be realized. The exchange of the tangential component of the velocity models the dynamic resistance to tension, and the exchange of the normal velocity characteristic is the dynamic resistance to bending. Together, this will accelerate the transmission of disturbances along the hair, which will improve dynamics, make the hair visually more connected and less springy.
In addition, there will also be a static tendency to straighten. Each point will seek to compensate for the fold of hair. If there is a bend at a point, a force proportional to the bend value and directed in the direction of decreasing the bend value will act on the point. Two points adjacent to the bend point will experience half the force in the opposite direction.
These interactions are enough to simulate the physics of hair, but we are not limited to it. It is necessary to add the interaction of the hair with solid objects. This makes practical sense. The point is not only that physical models as a rule include the interaction between different parallel-modeled entities, for example, liquids and solids. But also in the fact that in practical tasks, for example, in games, the GPU simulation must interact in real time with objects calculated on the CPU side. So I could not help but pay attention to such an interaction. Our hair will interact with solids, information about which will be transmitted in the video memory in each measure.
For simplicity, we will work only with round objects. On the CPU side, we will have several circle colliders from standard 2D unit physics. And the rule of interaction will be this: if the point of the hair is inside the solid, it is transferred outward, and the fraction directed towards the body is subtracted from the speed of such a point, and the same fraction is transferred to the body. We will not consider the absolute speed of the body, for simplicity.
These three points are too much connected to discuss them separately.
To describe the point from which many hairs are made, we use the following structure:
This structure is declared twice: on the CPU side and on the GPU side. For convenience. On the CPU side, we write the initial data, copy it to the GPU buffer, and then they are processed there. But it was possible to explain this structure only on the GPU side, if we do not need to transmit the initial data.
As for the parameters dummy1 and dummy2. In an article written by an engineer from nvidia, I read that it is better to keep the data of the video memory buffers in multiples of 128 bits. As this reduces the number of operations needed to calculate the offset.
The values of the remaining parameters, I think, are clear. Although, an attentive reader may ask: why is speed a type of float, and a change of speed is int? Short answer: because the change in speed is modified simultaneously by parallel threads, and to avoid errors in the calculations, you need to use a secure record. And the protected write function only works with integer variables. I will talk more about this below.
We have many points with which we model hair. Data on all points is stored in video memory and is accessible through the buffer interface:
In the shader code, we determine only its name and data type, and its size is set externally, from the side of the code executed on the processor.
How is the computer shader code structured, what is it all about? The code consists of kernels. This is the same as the methods, but each kernel runs in parallel on multiple cores. Therefore, for each, the number of flows is indicated in the form of a three-dimensional structure.
This is what an empty kernel looks like, in which there is no code, only the necessary external information:
The kernel has an input parameter id, which stores the three-dimensional index of the stream. This is very convenient, each thread knows its own index, which means that it can work with its own separate data unit.
On the processor side, the kernel is called like this:
These three digits “2, 2, 1” are connected with the line preceding the corresponding kernel:
These two triples of digits determine the number of threads, that is, the number of parallel kernel instances. You just need to multiply them: 8 * 4 * 1 * 2 * 2 * 1 = 128 threads.
Addressing flows will be on each axis. In this case, the x axis will be 8 * 2 = 16 units. On the y axis, 4 * 2 = 8 units. That is, if the kernel is called like this:
And on the shader side, the number of threads is set like this:
Then we will have (X * x) * (Y * y) * (Z * z) threads.
For example, suppose we need to process each pixel of a texture with a size of 256 x 256, and we want a separate stream to deal with each pixel. So, we can determine the number of threads as follows:
and on the shader side:
Inside the kernel, the id.x parameter will take values in the range [0, 255], the same as the id.y parameter
. So, here is the line:
each of 65536 pixels of the
id.xy texture will be painted white - this is the same as uint2 (id.x, id.y)
If this part related to the number of threads is unclear to someone, I advise you to go to the one I mentioned easier guide, and see how all this is used in practice to draw the Mandelbrot fractal using the simplest shader.
The shader text in the model we are considering contains several kernels, which in turn are launched on the CPU side in the Update () method. I will then review the text of each kernel, and first briefly explain what each of them does.
calc - the tangential and normal forces of interaction between the particles are calculated: the "spring" tension force pushes the particles along the line between them, and the "bending stiffness" forces the particles perpendicular to the line between adjacent particles; calculated forces are stored for each particle
velShare - particles exchange relative speeds. Tangential and fully inclusive - individually. Why highlight the tangential, if then there is still an exchange of full speed? The exchange of tangential velocity should be much more intense than normal, it should have a higher coefficient, so it had to be distinguished. Then why in the second case, I do not use the pure normal component, but use the full speed? To save on calculations. Changes in speed are recorded in the form of forces, As in the previous kernel.
interactionWithColliders - each point interacts with colliders, information about which is contained in the calcApply
buffer updated in each cycle - the forces calculated in previous kernels are added to the speed, and the speeds change the coordinates of the points
visInternodeLines - between the points, lines are drawn in a special buffer with a length of 1024 x 1024 (not yet on the texture)
pixelsToTexture - and here the values from the aforementioned are already converted to the colors of the pixels on the texture with the size [1024, 1024]
clearPixels - all the values of the intermediate buffer (in which we drew lines)
clearTexture is reset - oneThreadAction texture is cleared
- this kernel is executed in one single thread, it is needed to smoothly move the entire hair system to the place where we dragged it with the mouse. Smoothness is needed so that the system does not go abruptly from abrupt movement (as you remember, in our model the forces between particles are proportional to the square of the distance between them).
Now I will show how these kernels are launched from the side of the CPU code. But first, how to prepare the shader for launch.
Declare a variable:
We initialize it by specifying a file with shader text:
We set constants that are useful to us on the GPU side
We declare variables for the array that will store the data of the modeled points, and for the buffer, through whose interface we can read and write data to video memory
Initialize the buffer and write the array data to video memory.
For each kernel, set the used buffers so that the kernel can read and write data to this buffer.
When all the necessary buffers are created and installed for all shader kernels, you can run the kernels.
All kernels are launched from Update (). From FixedUpdate () they should not be launched (it will lag strongly), because the graphics pipeline is synchronized with Update ().
Kernels are launched in this sequence (I quote the entire code of the “doShaderStuff” method called in Update ()):
It immediately catches the eye that several kernels run 40 times per update. What for? So that with a small time step, the simulation works quickly in real time. And why should the time step be small? To reduce sampling error, that is, for system stability. And how and why does instability arise? If the step is large, and a large force acts on the point, then in one step the point flies away, the return force becomes even greater, and in the next step the point flies in the other direction even further. Result: the system goes peddling, all points fly back and forth with increasing amplitude. And with a small step, all the force and velocity curves are very smooth, because the errors greatly decrease with decreasing time step.
So instead of one big step, the system takes 40 small steps in each cycle, and thanks to this, it demonstrates high accuracy of calculations. Due to its high accuracy, it is possible to work with large forces of interaction without loss of stability. And great strength means that we have not sluggish, springy pasta in the model hanging out, trying to explode from a sudden movement, and durable hair spins cheerfully.
The data on the points with which we model the hair is stored in video memory in the form of a one-dimensional array, which we access through the buffer interface.
For the convenience of working with a one-dimensional buffer, we index the flows as follows: (x axis: number of hair * axis y: number of points in the hair). That is, we will have a two-dimensional array of streams, each of which will know its point by the stream index.
As you remember, the number of threads in which the kernel is executed is determined by the product of the parameters of the Dispatch () method and the parameters of the [numthreads ()] directive in the shader code.
In our case, all kernels that work with hair dots are preceded by the [numthreads (16.8.1)] directive. Therefore, the parameters of the Dispatch () method must be such that the product gives the number of threads no less than we need to process the entire array of points. In the code, we calculate the x and y parameters of the Dispatch () method:
The relationship between the parameters [numthreads ()] and Dispatch () stems from the architecture of graphic computers. The first is the number of threads in the group. The second is the number of thread groups. Their ratio affects the speed of work. If we need 1024 streams along the x axis, it is better to make 32 groups of 32 streams than 1 group of 1024 streams. Why? To answer this question, you need to talk a lot about the architecture of the GPU, let us leave this too deep topic untouched.
So, 40 times in an update, we launch in turn kernels that calculate the change in the speed of points and change their speed and coordinates. Let's look at the code for each kernel. Everything is quite simple there, you just need to learn a couple of specific features.
Kernel calc calculates the change in speed of the points. The points in the hairNodesBuffer buffer are arranged in turn, first the first point of the first hair, then the second, and so on to the last. Then immediately the first point of the second hair, and so on through all the hair, to the end of the buffer. We remember that the kernel has an id parameter, and in our case id.x indicates the hair number, and id.y - the point number. And here is how we get access to the data points:
Here, the value nNodesPerHair is the constant that we set on the CPU side when initializing the shader. The data from the buffer is copied to the local variables node and node2 because accessing the buffer data may require more kernel cycles than accessing the local variable. The algorithm itself is as follows: for each point, if it is not the last in the hair, we calculate the force acting between it and the next point. Based on this force, we record a change in speed at each of the points.
Here is an important feature of parallel computing: each stream modifies two points, the current and the next, which means that two parallel flows modify each point. Unprotected writing to variables common to parallel threads is fraught with data loss. If you use the usual increment:
then recording can happen simultaneously, like this: the first stream will copy the original value, add one to it, but before it writes the value back to the memory cell, the second stream will take the original value. Then the first thread will write the value increased by one back. After which the second stream will add its unit and write the increased value back. Result: although two threads were added one at a time, the variable increased by only one unit. To avoid this situation, use a secure recording. HLSL has several functions for securely modifying generic variables. They guarantee that the data will not be lost and the contribution of each stream will be taken into account.
A small problem is that these functions only work with integer variables. And that is why in the structure describing the state of the point we use the dvx and dvy parameters of type int. To be able to write to them using protected functions and not lose data. But in order not to lose accuracy on rounding, we determined the factors in advance. One translates the float to int, the other back. So we use the entire range of int-values, and we do not lose in accuracy (we lose, of course, but negligible).
A protected record looks like this:
Here F_TO_I is the mentioned coefficient for the projection of float on int, dv is the force vector of the influence of the second particle on the first through the spring connection. And dvFlex is a rectifying force. "(int)" needs to be added because InterlockedAdd () is overloaded for int and uint types, and float is interpreted as uint by default.
The velShare Kernel is similar to the previous one, it also modifies the dvx and dvy parameters of two adjacent points, but instead of calculating the forces, the diffusion of the relative speed is calculated.
In the interactionWithColliders kernel, points do not interact with each other, here each point runs through all the colliders of the solid buffer (which we update in each update). That is, each stream writes only to one particle, there is no danger of simultaneous recording, and therefore, instead of InterlockedAdd (), we can directly change the speed of the particle. But at the same time, our model implies that the points transmit momentum to the collider. This means that parallel streams can simultaneously change the momentum of the same collider, which means we use a protected recording option.
Only here you need to understand: when we project a float on an int, the integer and fractional parts compete with us. Accuracy competes with a range of magnitude. For the case of interaction of points, we chose a coefficient that admits a sufficient dispersion of magnitude for us, and the rest was allowed for accuracy. But this coefficient is not suitable for transferring momentum to the collider, because at the same time hundreds of points can add their momentum in one direction, and therefore it is necessary to sacrifice accuracy in favor of the ability to accommodate a large number. So with a protected record, we do not use the coefficient F_TO_I, but use a smaller coefficient.
After all the interactions of the points are calculated, in the calcApply kernel we add momentum to speed, and speed to coordinates. In addition, in this kernel, each root (first) point of the hair is fixed in a certain place relative to the current position of the entire hair system. Even in this kernel, the contribution of gravity is added to the vertical velocity component. Plus, “braking” about air is realized, that is, the absolute value of the speed of each point is multiplied by a coefficient slightly less than unity.
Note that in the calcApply kernel, speed affects the coordinates through the dPosRate coefficient. It determines the size of the simulation step. This coefficient is set on the CPU side and is stored in a variable, which I called "simulationSpeed". The larger this parameter, the faster the system will evolve over time. But the lower will be the accuracy of the calculation. The accuracy of the calculation, I repeat, limits the magnitude of the forces, since for large forces and low accuracy, the magnitude of the error is so great that it determines the behavior of the model. We took the simulation speed rather low, this gives us greater accuracy, so we can afford greater forces, which means more realistic model behavior.
For the magnitude of the forces, there is a coefficient that relates the effect of the pulse on the speed - dVelRate. This coefficient is large, it is set on the CPU side and is called "strengthOfForces".
I repeat that in all the mentioned kernels the number of threads is equal to the number of points, one thread is responsible for processing one point. And this is a good practice. We pay nothing for the number of threads, there can be any number of them (in shader model 5.0 - no more than 1024 along the x and y axes and no more than 64 along the z axis). In the tradition of parallel computing, it is better to avoid using loops to perform a single operation in one thread in relation to several units of data, it is better to make as many threads as are necessary to implement the principle of “one unit of data - one stream”.
We return to the doShaderStuff () method on the side of the CPU code. After completing a cycle of 40 steps of calculating the hair model, we read the data of the collider:
You may recall that on the GPU side, pulses from the hair side are recorded in the buffer with collider data, and we use them on the CPU side to apply force to the rigidbody. Note that the force on the rigidbody is applied in the FixedUpdate () method, since it is synchronized with physics. In this case, the pulse data is updated in Update (). So, under the influence of various factors, several FixedUpdate () and vice versa can happen in one Update (). That is, there is no absolute accuracy in the effect of hair on the collider, some of the data can be overwritten before it is affected, and other data can be affected twice. You can take measures to prevent this from happening, but these measures are not taken in the program under consideration.
It is also worth noting that the GetData () method pauses the graphics pipeline, which causes a noticeable slowdown. Unfortunately, the asynchronous version of this method in the unit has not yet been implemented, although it is rumored that in 2018 it will appear. In the meantime, you need to understand that if in your task you need to copy data from the GPU to the CPU, the program will work 20-30% slower. At the same time, the SetData () method does not have such an effect, it works quickly.
The remaining kernels launched in the doShaderStuff () method are associated only with the visualization of the hair system.
Consider everything related to visualization.
On the CPU side, we declare a RenderTexture variable, do not forget to set enableRandomWrite = true, and use it as mainTexture in the material of the Image UI component.
And then for each kernel that should write to this texture, we call the SetTexture () method to bind our RenderTexture object to a variable on the shader side:
On the shader side, we have declared a variable of the type RWTexture2D, through which we set the texture pixel colors:
Now consider the texture cleaning kernel that is called before writing color pixels to it:
This kernel runs like this:
We see that we have 1024 x 1024 streams, per pixel per stream. Which is convenient: just use the id.xy parameter to address the pixel.
How exactly is hair drawn? I decided to make the hair translucent, so that when they intersect, the color is more saturated, which implies the need to use a secure recording, since two lines can be drawn at the same pixel at the same time, since, like in the kernels already considered, all points will be executed simultaneously in the number of threads equal to the number of points. Drawing itself is trivial: from each point we draw a line to the next point. There are special algorithms for choosing the set of square pixels swept by the line, but I decided to go the simple way: the line is drawn by moving in small steps along the line between two points.
Since increment is used, I write color data to a buffer, not a texture. The texture for some reason is not readable, although it seems to be.
After the visInternodeLines kernel has drawn all the lines, we copy the pixels from the buffer to the texture. I did not use any colors, only gradations of gray are drawn. If I needed color, instead of the RWStructuredBuffer buffer, I would use RWStructuredBuffer or I could write 4 color parameters to one uint.
By the way, this method with RenderTexture does not work on poppies, and I could not get an answer to the question “why” on the forum.
There are other methods for visualizing data from the compute shader, but I must confess that I have not studied them yet.
After the “pixelsToTexture” kernel modified the texture, images of flying hair appear on our screen.
I talked about all the pieces of code regarding GPU computing. There is a lot of information in the manual, and it can be difficult to grasp it at once. If you plan to experiment in this area, I advise you to write a simple program from scratch in order to consolidate knowledge through practice. Think of it as homework. Performing it will be simple and useful.
Write a shader with one kernel squaring all the numbers in a large array. On the CPU side, prepare an array, write it to the shader buffer, start the kernel, then get the information from the video memory and check if the numbers are squared.
Here is a project for Unity3D, on the explanation of which the guide is built. You need to download and open it in Unity: a
link to the Unity project.
Who will understand this guide? Those who use Unity3D or at least know C # or C ++. The shader is written in HLSL, a close syntactic cousin of C ++.
Who will benefit from this guide? Experienced programmers who want to learn how to use the GPU for computing. But even an inexperienced but diligent programmer will easily understand everything.
Why use a graphics card for computing? For parallel tasks, its performance is 10-100 times higher than that of the processor. That is, everyone in the computer has a small supercomputer with a convenient API, it makes sense to use it in suitable cases.
Is this huge performance really needed? Yes, processor speed is often a limiting factor. For example, when you need to perform the same operations on large data sets. But it is precisely such tasks that can be easily parallelized. In addition, developers often refuse decisions because of their computing capacity, and entire areas in the algorithm space remain unexplored. For example, you can do the coolest physics in games, if you properly load the GPU.
And what, with a video card you can now simply solve problems with brute force? The demand for optimization does not depend on the performance of iron. There is no such supercomputer that could not be tightly loaded with inefficient code.
Why exactly compute shader? Why not opencl or cuda? Cuda only works on nvidia hardware, and I don’t know opencl. Unity can build in any API, including opengl core. On Macs and on an Android computer, shaders work, on Linux, it seems, too (although I have not tried it). Although, each API has limitations that should be considered. For example, on Metal you cannot make more than 256 threads along one axis (In DX10 - 1024). And the android API will not be able to use more than 4 buffers per kernel (in DX10 - 8, in DX11 - even more).
Why physical simulation? This is a computationally intensive task, and is well suited for parallel computing. In addition, the task is in demand. Game devs can implement interesting physics in games, students can create experimental models for term papers, engineers and scientists can do calculations on models.
And why exactly the hair model? I wanted to take a simple puzzle, but at the same time covering the main problems.
How to use this manual? It is best to download the source code, open it and read it as you progress through the manual. I will explain in detail all the main lines, although I will not explain absolutely every line, the meaning of most of them is obvious. There are no complex algorithms in the text, there is only the use of the interface of classes that serve computing on the GPU. And on the side of the shader code, there is nothing but reading the data, performing simple mathematical operations on them, and recording the results. But if something is not clear - be sure to ask, I will answer everything in kamenty.
And now for those who have absolutely no idea about using compute shaders, I suggest taking a step to the side and moving on to a very simple guide, which is devoted to the basics of using computer shaders. I advise you to start with it in order to better understand the essence and get used to the practice of GPU computing with an extremely simple example. And then come back here and continue. And those who are at least somehow familiar with computer shaders, let them boldly read on.
If you want to make a physical model calculated on a GPU from scratch, then this problem can be divided into 4 parts:
- a mathematical model of the phenomenon
- an algorithm for parallel model calculation
- a shader code
- preparing and running a shader in a unit
Mathematical model
The strengths of video cards are that they can apply one operation simultaneously to many objects. Therefore, a hair model can be made as a chain of dots, each of which interacts with two neighbors. The interaction between the points is based on the spring principle: k * (S0-S) ^ n, where S0 is the equilibrium distance, S is the current distance. In reality, the hair does not look like a spring, it is perceived as inextensible. This means that the spring in the model must be made sufficiently rigid. It is better to increase the stiffness of the spring by increasing n, because the degree increases the curvature of the curve in the vicinity of the equilibrium, which reduces the backlash and reduces the effect of “rubberiness” of the hair. I took n = 2, and we will talk about the value of the coefficient k below.
In addition to the elastic force between the points, diffusion of relative velocities or one-dimensional viscosity will be realized. The exchange of the tangential component of the velocity models the dynamic resistance to tension, and the exchange of the normal velocity characteristic is the dynamic resistance to bending. Together, this will accelerate the transmission of disturbances along the hair, which will improve dynamics, make the hair visually more connected and less springy.
In addition, there will also be a static tendency to straighten. Each point will seek to compensate for the fold of hair. If there is a bend at a point, a force proportional to the bend value and directed in the direction of decreasing the bend value will act on the point. Two points adjacent to the bend point will experience half the force in the opposite direction.
These interactions are enough to simulate the physics of hair, but we are not limited to it. It is necessary to add the interaction of the hair with solid objects. This makes practical sense. The point is not only that physical models as a rule include the interaction between different parallel-modeled entities, for example, liquids and solids. But also in the fact that in practical tasks, for example, in games, the GPU simulation must interact in real time with objects calculated on the CPU side. So I could not help but pay attention to such an interaction. Our hair will interact with solids, information about which will be transmitted in the video memory in each measure.
For simplicity, we will work only with round objects. On the CPU side, we will have several circle colliders from standard 2D unit physics. And the rule of interaction will be this: if the point of the hair is inside the solid, it is transferred outward, and the fraction directed towards the body is subtracted from the speed of such a point, and the same fraction is transferred to the body. We will not consider the absolute speed of the body, for simplicity.
Algorithm, code and shader preparation
These three points are too much connected to discuss them separately.
To describe the point from which many hairs are made, we use the following structure:
struct hairNode{
float x; // позиция в двумерном пространстве
float y; //
float vx; // скорость
float vy; //
int dvx; // сумма сил - изменение скорости на данном шаге
int dvy; //
int dummy1; // эти два параметра добавлены для кратности 128 битам
int dummy2; //
}This structure is declared twice: on the CPU side and on the GPU side. For convenience. On the CPU side, we write the initial data, copy it to the GPU buffer, and then they are processed there. But it was possible to explain this structure only on the GPU side, if we do not need to transmit the initial data.
As for the parameters dummy1 and dummy2. In an article written by an engineer from nvidia, I read that it is better to keep the data of the video memory buffers in multiples of 128 bits. As this reduces the number of operations needed to calculate the offset.
The values of the remaining parameters, I think, are clear. Although, an attentive reader may ask: why is speed a type of float, and a change of speed is int? Short answer: because the change in speed is modified simultaneously by parallel threads, and to avoid errors in the calculations, you need to use a secure record. And the protected write function only works with integer variables. I will talk more about this below.
We have many points with which we model hair. Data on all points is stored in video memory and is accessible through the buffer interface:
RWStructuredBuffer hairNodesBuffer; In the shader code, we determine only its name and data type, and its size is set externally, from the side of the code executed on the processor.
How is the computer shader code structured, what is it all about? The code consists of kernels. This is the same as the methods, but each kernel runs in parallel on multiple cores. Therefore, for each, the number of flows is indicated in the form of a three-dimensional structure.
This is what an empty kernel looks like, in which there is no code, only the necessary external information:
#pragma kernel kernelName
[numthreads(8,4,1)]
void kernelName (uint3 id : SV_DispatchThreadID){
// здесь должен быть код
}The kernel has an input parameter id, which stores the three-dimensional index of the stream. This is very convenient, each thread knows its own index, which means that it can work with its own separate data unit.
On the processor side, the kernel is called like this:
shaderInstance.Dispatch(kernelIndex, 2, 2, 1);These three digits “2, 2, 1” are connected with the line preceding the corresponding kernel:
[numthreads(8,4,1)]These two triples of digits determine the number of threads, that is, the number of parallel kernel instances. You just need to multiply them: 8 * 4 * 1 * 2 * 2 * 1 = 128 threads.
Addressing flows will be on each axis. In this case, the x axis will be 8 * 2 = 16 units. On the y axis, 4 * 2 = 8 units. That is, if the kernel is called like this:
ComputeShader.Dispatch(kernelIndex, X, Y, Z);And on the shader side, the number of threads is set like this:
[numthreads(x,y,z)]Then we will have (X * x) * (Y * y) * (Z * z) threads.
For example, suppose we need to process each pixel of a texture with a size of 256 x 256, and we want a separate stream to deal with each pixel. So, we can determine the number of threads as follows:
Dispatch(kernelIndex, 16, 16, 1);and on the shader side:
[numthreads(16,16,1)]Inside the kernel, the id.x parameter will take values in the range [0, 255], the same as the id.y parameter
. So, here is the line:
texture[id.xy]=float4(1, 1, 1, 1);each of 65536 pixels of the
id.xy texture will be painted white - this is the same as uint2 (id.x, id.y)
If this part related to the number of threads is unclear to someone, I advise you to go to the one I mentioned easier guide, and see how all this is used in practice to draw the Mandelbrot fractal using the simplest shader.
The shader text in the model we are considering contains several kernels, which in turn are launched on the CPU side in the Update () method. I will then review the text of each kernel, and first briefly explain what each of them does.
calc - the tangential and normal forces of interaction between the particles are calculated: the "spring" tension force pushes the particles along the line between them, and the "bending stiffness" forces the particles perpendicular to the line between adjacent particles; calculated forces are stored for each particle
velShare - particles exchange relative speeds. Tangential and fully inclusive - individually. Why highlight the tangential, if then there is still an exchange of full speed? The exchange of tangential velocity should be much more intense than normal, it should have a higher coefficient, so it had to be distinguished. Then why in the second case, I do not use the pure normal component, but use the full speed? To save on calculations. Changes in speed are recorded in the form of forces, As in the previous kernel.
interactionWithColliders - each point interacts with colliders, information about which is contained in the calcApply
buffer updated in each cycle - the forces calculated in previous kernels are added to the speed, and the speeds change the coordinates of the points
visInternodeLines - between the points, lines are drawn in a special buffer with a length of 1024 x 1024 (not yet on the texture)
pixelsToTexture - and here the values from the aforementioned are already converted to the colors of the pixels on the texture with the size [1024, 1024]
clearPixels - all the values of the intermediate buffer (in which we drew lines)
clearTexture is reset - oneThreadAction texture is cleared
- this kernel is executed in one single thread, it is needed to smoothly move the entire hair system to the place where we dragged it with the mouse. Smoothness is needed so that the system does not go abruptly from abrupt movement (as you remember, in our model the forces between particles are proportional to the square of the distance between them).
CPU side
Now I will show how these kernels are launched from the side of the CPU code. But first, how to prepare the shader for launch.
Declare a variable:
ComputeShader _shader;We initialize it by specifying a file with shader text:
_shader = Resources.Load("shader"); We set constants that are useful to us on the GPU side
// переменные nodesPerHair и nHairs уже инициализированы
_shader.SetInt("nNodsPerHair", nodesPerHair);
_shader.SetInt("nHairs", nHairs);We declare variables for the array that will store the data of the modeled points, and for the buffer, through whose interface we can read and write data to video memory
hairNode[] hairNodesArray;
ComputeBuffer hairNodesBuffer;Initialize the buffer and write the array data to video memory.
// hairNodesArray уже инициализирован
hairNodesBuffer = new ComputeBuffer(hairNodesArray.Length, 4 * 8);
hairNodesBuffer.SetData(hairNodesArray);For each kernel, set the used buffers so that the kernel can read and write data to this buffer.
kiCalc = _shader.FindKernel("calc");
_shader.SetBuffer(kiCalc, "hairNodesBuffer", hairNodesBuffer);When all the necessary buffers are created and installed for all shader kernels, you can run the kernels.
All kernels are launched from Update (). From FixedUpdate () they should not be launched (it will lag strongly), because the graphics pipeline is synchronized with Update ().
Kernels are launched in this sequence (I quote the entire code of the “doShaderStuff” method called in Update ()):
void doShaderStuff(){
int i, nHairThreadGroups, nNodesThreadGroups;
nHairThreadGroups = (nHairs - 1) / 16 + 1;
nNodesThreadGroups = (nodesPerHair - 1) / 8 + 1;
_shader.SetFloats("pivotDestination", pivotPosition);
circleCollidersBuffer.SetData(circleCollidersArray);
i = 0;
while (i < 40) {
_shader.Dispatch(kiVelShare, nHairThreadGroups, nNodesThreadGroups, 1);
_shader.Dispatch(kiCalc, nHairThreadGroups, nNodesThreadGroups, 1);
_shader.Dispatch(kiInteractionWithColliders, nHairThreadGroups, nNodesThreadGroups, 1);
_shader.Dispatch(kiCalcApply, nHairThreadGroups, nNodesThreadGroups, 1);
_shader.Dispatch(kiOneThreadAction, 1, 1, 1);
i++;
}
circleCollidersBuffer.GetData(circleCollidersArray);
_shader.Dispatch(kiVisInternodeLines, nHairThreadGroups, nNodesThreadGroups, 1);
_shader.Dispatch(kiClearTexture, 32, 32, 1);
_shader.Dispatch(kiPixelsToTexture, 32, 32, 1);
_shader.Dispatch(kiClearPixels, 32, 32, 1);
}
It immediately catches the eye that several kernels run 40 times per update. What for? So that with a small time step, the simulation works quickly in real time. And why should the time step be small? To reduce sampling error, that is, for system stability. And how and why does instability arise? If the step is large, and a large force acts on the point, then in one step the point flies away, the return force becomes even greater, and in the next step the point flies in the other direction even further. Result: the system goes peddling, all points fly back and forth with increasing amplitude. And with a small step, all the force and velocity curves are very smooth, because the errors greatly decrease with decreasing time step.
So instead of one big step, the system takes 40 small steps in each cycle, and thanks to this, it demonstrates high accuracy of calculations. Due to its high accuracy, it is possible to work with large forces of interaction without loss of stability. And great strength means that we have not sluggish, springy pasta in the model hanging out, trying to explode from a sudden movement, and durable hair spins cheerfully.
The data on the points with which we model the hair is stored in video memory in the form of a one-dimensional array, which we access through the buffer interface.
For the convenience of working with a one-dimensional buffer, we index the flows as follows: (x axis: number of hair * axis y: number of points in the hair). That is, we will have a two-dimensional array of streams, each of which will know its point by the stream index.
As you remember, the number of threads in which the kernel is executed is determined by the product of the parameters of the Dispatch () method and the parameters of the [numthreads ()] directive in the shader code.
In our case, all kernels that work with hair dots are preceded by the [numthreads (16.8.1)] directive. Therefore, the parameters of the Dispatch () method must be such that the product gives the number of threads no less than we need to process the entire array of points. In the code, we calculate the x and y parameters of the Dispatch () method:
nHairThreadGroups = (nHairs - 1) / 16 + 1;
nNodesThreadGroups = (nodesPerHair - 1) / 8 + 1;The relationship between the parameters [numthreads ()] and Dispatch () stems from the architecture of graphic computers. The first is the number of threads in the group. The second is the number of thread groups. Their ratio affects the speed of work. If we need 1024 streams along the x axis, it is better to make 32 groups of 32 streams than 1 group of 1024 streams. Why? To answer this question, you need to talk a lot about the architecture of the GPU, let us leave this too deep topic untouched.
GPU Details
So, 40 times in an update, we launch in turn kernels that calculate the change in the speed of points and change their speed and coordinates. Let's look at the code for each kernel. Everything is quite simple there, you just need to learn a couple of specific features.
Kernel calc calculates the change in speed of the points. The points in the hairNodesBuffer buffer are arranged in turn, first the first point of the first hair, then the second, and so on to the last. Then immediately the first point of the second hair, and so on through all the hair, to the end of the buffer. We remember that the kernel has an id parameter, and in our case id.x indicates the hair number, and id.y - the point number. And here is how we get access to the data points:
int nodeIndex, nodeIndex2;
hairNode node, node2;
nodeIndex = id.x * nNodesPerHair + id.y;
nodeIndex2 = nodeIndex + 1;
node = hairNodesBuffer[nodeIndex];
node2 = hairNodesBuffer[nodeIndex2];Here, the value nNodesPerHair is the constant that we set on the CPU side when initializing the shader. The data from the buffer is copied to the local variables node and node2 because accessing the buffer data may require more kernel cycles than accessing the local variable. The algorithm itself is as follows: for each point, if it is not the last in the hair, we calculate the force acting between it and the next point. Based on this force, we record a change in speed at each of the points.
Here is an important feature of parallel computing: each stream modifies two points, the current and the next, which means that two parallel flows modify each point. Unprotected writing to variables common to parallel threads is fraught with data loss. If you use the usual increment:
variable += value;then recording can happen simultaneously, like this: the first stream will copy the original value, add one to it, but before it writes the value back to the memory cell, the second stream will take the original value. Then the first thread will write the value increased by one back. After which the second stream will add its unit and write the increased value back. Result: although two threads were added one at a time, the variable increased by only one unit. To avoid this situation, use a secure recording. HLSL has several functions for securely modifying generic variables. They guarantee that the data will not be lost and the contribution of each stream will be taken into account.
A small problem is that these functions only work with integer variables. And that is why in the structure describing the state of the point we use the dvx and dvy parameters of type int. To be able to write to them using protected functions and not lose data. But in order not to lose accuracy on rounding, we determined the factors in advance. One translates the float to int, the other back. So we use the entire range of int-values, and we do not lose in accuracy (we lose, of course, but negligible).
A protected record looks like this:
InterlockedAdd(hairNodesBuffer[nodeIndex].dvx, (int)(F_TO_I * (dv.x + 2 * dvFlex.x)));
InterlockedAdd(hairNodesBuffer[nodeIndex].dvy, (int)(F_TO_I * (dv.y + 2 * dvFlex.y)));
InterlockedAdd(hairNodesBuffer[nodeIndex2].dvx, (int)(F_TO_I * (-dv.x - dvFlex.x)));
InterlockedAdd(hairNodesBuffer[nodeIndex2].dvy, (int)(F_TO_I * (-dv.y - dvFlex.y)));Here F_TO_I is the mentioned coefficient for the projection of float on int, dv is the force vector of the influence of the second particle on the first through the spring connection. And dvFlex is a rectifying force. "(int)" needs to be added because InterlockedAdd () is overloaded for int and uint types, and float is interpreted as uint by default.
The velShare Kernel is similar to the previous one, it also modifies the dvx and dvy parameters of two adjacent points, but instead of calculating the forces, the diffusion of the relative speed is calculated.
In the interactionWithColliders kernel, points do not interact with each other, here each point runs through all the colliders of the solid buffer (which we update in each update). That is, each stream writes only to one particle, there is no danger of simultaneous recording, and therefore, instead of InterlockedAdd (), we can directly change the speed of the particle. But at the same time, our model implies that the points transmit momentum to the collider. This means that parallel streams can simultaneously change the momentum of the same collider, which means we use a protected recording option.
Only here you need to understand: when we project a float on an int, the integer and fractional parts compete with us. Accuracy competes with a range of magnitude. For the case of interaction of points, we chose a coefficient that admits a sufficient dispersion of magnitude for us, and the rest was allowed for accuracy. But this coefficient is not suitable for transferring momentum to the collider, because at the same time hundreds of points can add their momentum in one direction, and therefore it is necessary to sacrifice accuracy in favor of the ability to accommodate a large number. So with a protected record, we do not use the coefficient F_TO_I, but use a smaller coefficient.
After all the interactions of the points are calculated, in the calcApply kernel we add momentum to speed, and speed to coordinates. In addition, in this kernel, each root (first) point of the hair is fixed in a certain place relative to the current position of the entire hair system. Even in this kernel, the contribution of gravity is added to the vertical velocity component. Plus, “braking” about air is realized, that is, the absolute value of the speed of each point is multiplied by a coefficient slightly less than unity.
Note that in the calcApply kernel, speed affects the coordinates through the dPosRate coefficient. It determines the size of the simulation step. This coefficient is set on the CPU side and is stored in a variable, which I called "simulationSpeed". The larger this parameter, the faster the system will evolve over time. But the lower will be the accuracy of the calculation. The accuracy of the calculation, I repeat, limits the magnitude of the forces, since for large forces and low accuracy, the magnitude of the error is so great that it determines the behavior of the model. We took the simulation speed rather low, this gives us greater accuracy, so we can afford greater forces, which means more realistic model behavior.
For the magnitude of the forces, there is a coefficient that relates the effect of the pulse on the speed - dVelRate. This coefficient is large, it is set on the CPU side and is called "strengthOfForces".
I repeat that in all the mentioned kernels the number of threads is equal to the number of points, one thread is responsible for processing one point. And this is a good practice. We pay nothing for the number of threads, there can be any number of them (in shader model 5.0 - no more than 1024 along the x and y axes and no more than 64 along the z axis). In the tradition of parallel computing, it is better to avoid using loops to perform a single operation in one thread in relation to several units of data, it is better to make as many threads as are necessary to implement the principle of “one unit of data - one stream”.
We return to the doShaderStuff () method on the side of the CPU code. After completing a cycle of 40 steps of calculating the hair model, we read the data of the collider:
circleCollidersBuffer.GetData(circleCollidersArray);You may recall that on the GPU side, pulses from the hair side are recorded in the buffer with collider data, and we use them on the CPU side to apply force to the rigidbody. Note that the force on the rigidbody is applied in the FixedUpdate () method, since it is synchronized with physics. In this case, the pulse data is updated in Update (). So, under the influence of various factors, several FixedUpdate () and vice versa can happen in one Update (). That is, there is no absolute accuracy in the effect of hair on the collider, some of the data can be overwritten before it is affected, and other data can be affected twice. You can take measures to prevent this from happening, but these measures are not taken in the program under consideration.
It is also worth noting that the GetData () method pauses the graphics pipeline, which causes a noticeable slowdown. Unfortunately, the asynchronous version of this method in the unit has not yet been implemented, although it is rumored that in 2018 it will appear. In the meantime, you need to understand that if in your task you need to copy data from the GPU to the CPU, the program will work 20-30% slower. At the same time, the SetData () method does not have such an effect, it works quickly.
Visualization
The remaining kernels launched in the doShaderStuff () method are associated only with the visualization of the hair system.
Consider everything related to visualization.
On the CPU side, we declare a RenderTexture variable, do not forget to set enableRandomWrite = true, and use it as mainTexture in the material of the Image UI component.
And then for each kernel that should write to this texture, we call the SetTexture () method to bind our RenderTexture object to a variable on the shader side:
RenderTexture renderTexture;
renderTexture = new RenderTexture(1024, 1024, 32);
renderTexture.enableRandomWrite = true;
renderTexture.Create();
GameObject.Find("canvas/image").GetComponent().material.mainTexture = renderTexture;
_shader.SetTexture(kiPixelsToTexture, "renderTexture", renderTexture); On the shader side, we have declared a variable of the type RWTexture2D, through which we set the texture pixel colors:
RWTexture2D renderTexture; Now consider the texture cleaning kernel that is called before writing color pixels to it:
#pragma kernel clearTexture
[numthreads(32,32,1)]
void clearTexture (uint3 id : SV_DispatchThreadID){
renderTexture[id.xy] = float4(0, 0, 0, 0);
}This kernel runs like this:
_shader.Dispatch(kiClearTexture, 32, 32, 1);We see that we have 1024 x 1024 streams, per pixel per stream. Which is convenient: just use the id.xy parameter to address the pixel.
How exactly is hair drawn? I decided to make the hair translucent, so that when they intersect, the color is more saturated, which implies the need to use a secure recording, since two lines can be drawn at the same pixel at the same time, since, like in the kernels already considered, all points will be executed simultaneously in the number of threads equal to the number of points. Drawing itself is trivial: from each point we draw a line to the next point. There are special algorithms for choosing the set of square pixels swept by the line, but I decided to go the simple way: the line is drawn by moving in small steps along the line between two points.
Since increment is used, I write color data to a buffer, not a texture. The texture for some reason is not readable, although it seems to be.
After the visInternodeLines kernel has drawn all the lines, we copy the pixels from the buffer to the texture. I did not use any colors, only gradations of gray are drawn. If I needed color, instead of the RWStructuredBuffer buffer, I would use RWStructuredBuffer or I could write 4 color parameters to one uint.
By the way, this method with RenderTexture does not work on poppies, and I could not get an answer to the question “why” on the forum.
There are other methods for visualizing data from the compute shader, but I must confess that I have not studied them yet.
After the “pixelsToTexture” kernel modified the texture, images of flying hair appear on our screen.
I talked about all the pieces of code regarding GPU computing. There is a lot of information in the manual, and it can be difficult to grasp it at once. If you plan to experiment in this area, I advise you to write a simple program from scratch in order to consolidate knowledge through practice. Think of it as homework. Performing it will be simple and useful.
Write a shader with one kernel squaring all the numbers in a large array. On the CPU side, prepare an array, write it to the shader buffer, start the kernel, then get the information from the video memory and check if the numbers are squared.
