Approximating rigidbody motion in a viscous fluid by scripting additional forces

Hello, I'm relatively new to Unity scripting and trying to do some custom scripting working with Unity's PhysX physics system.

Basically, I want to simulate rigidbody cube primitives with an additional drag force to approximate an underwater environment, similar to Karl Sims' evolved virtual creatures swimming (first 30 secs):

From his 1994 paper:
"A viscosity effect is used for the simulations in underwater environments. For each exposed moving surface, a viscous force resists the normal component of its velocity, proportional to its surface area and normal velocity magnitude. This is a simple approximation that does not include the motion of the fluid itself, but is still sufficient for simulating realistic looking swimming and paddling dynamics. "

So in Unity, for each rigidbody, I want to iterate through each face, measure its velocity and surface area, and apply a force proportional to that in the negative normal direction. Sounds simple until I try to actually implement it.
--------------------------------------
CODE:

Code (CSharp):

1. using UnityEngine;
2. using System.Collections;
3.  
4. public class sub_box_fluid_drag_script : MonoBehaviour {
5.     // Debug shapes! (using them as a visual representation of where/how much force is applied)
6.     public GameObject db_zpos;
7.     public GameObject db_ypos;
8.     public GameObject db_xpos;
9.  
10.     public float viscosityDrag = 1f;
11.  
12.     private Rigidbody rigidBod;
13.  
14.     // Surface Areas for each pair of faces (neg x will be same as pos x):
15.     private float sa_x;
16.     private float sa_y;
17.     private float sa_z;
18.  
19.     // Use this for initialization
20.     void Start () {
21.         rigidBod = GetComponent<Rigidbody> ();
22.  
23.         // Calculate surface areas for each face:
24.         sa_x = transform.localScale.y * transform.localScale.z;
25.         sa_y = transform.localScale.x * transform.localScale.z;
26.         sa_z = transform.localScale.x * transform.localScale.y;
27.     }
28.    
29.     // Update is called once per frame
30.     void Update () {
31.    
32.     }
33.  
34.     void FixedUpdate () {
35.         // Find centers of each of box's faces
36.         Vector3 xpos_face_center = (transform.right * transform.localScale.x / 2) + transform.position;
37.         Vector3 ypos_face_center = (transform.up * transform.localScale.y / 2) + transform.position;
38.         Vector3 zpos_face_center = (transform.forward * transform.localScale.z / 2) + transform.position;
39.         Vector3 xneg_face_center = -(transform.right * transform.localScale.x / 2) + transform.position;
40.         Vector3 yneg_face_center = -(transform.up * transform.localScale.y / 2) + transform.position;
41.         Vector3 zneg_face_center = -(transform.forward * transform.localScale.z / 2) + transform.position;
42.  
43.         //=== FOR EACH FACE of rigidbody box: ----------------------------------------
44.         //=== Get Velocity: ---------------------------------------------
45.         //=== Apply Opposing Force ----------------------------------------
46.  
47.         // FRONT (posZ):
48.         Vector3 pointVelPosZ = rigidBod.GetPointVelocity (zpos_face_center); // Get velocity of face's center (doesn't catch torque around center of mass)
49.         Vector3 fluidDragVecPosZ = -transform.forward.normalized *     // in the direction opposite the face's normal
50.                                     Vector3.Dot (transform.forward, pointVelPosZ)   // get the proportion of the velocity vector in the direction of face's normal
51.                                     * sa_z * viscosityDrag;   // multiplied by face's surface area, and user-defined multiplier
52.         rigidBod.AddForceAtPosition (fluidDragVecPosZ, zpos_face_center);  // Apply force at face's center, in the direction opposite the face normal
53.         // BACK (negZ):
54.         Vector3 pointVelNegZ = rigidBod.GetPointVelocity (zneg_face_center);
55.         Vector3 fluidDragVecNegZ = transform.forward.normalized * Vector3.Dot (-transform.forward, pointVelNegZ) * sa_z * viscosityDrag;
56.         rigidBod.AddForceAtPosition (fluidDragVecNegZ, zneg_face_center);
57.  
58.         // TOP (posY):
59.         Vector3 pointVelPosY = rigidBod.GetPointVelocity (ypos_face_center);
60.         Vector3 fluidDragVecPosY = -transform.up.normalized * Vector3.Dot (transform.up, pointVelPosY) * sa_y * viscosityDrag;
61.         rigidBod.AddForceAtPosition (fluidDragVecPosY, ypos_face_center);
62.         // BOTTOM (negY):
63.         Vector3 pointVelNegY = rigidBod.GetPointVelocity (yneg_face_center);
64.         Vector3 fluidDragVecNegY = transform.up.normalized * Vector3.Dot (-transform.up, pointVelNegY) * sa_y * viscosityDrag;
65.         rigidBod.AddForceAtPosition (fluidDragVecNegY, yneg_face_center);
66.  
67.         // RIGHT (posX):
68.         Vector3 pointVelPosX = rigidBod.GetPointVelocity (xpos_face_center);
69.         Vector3 fluidDragVecPosX = -transform.right.normalized * Vector3.Dot (transform.right, pointVelPosX) * sa_x * viscosityDrag;
70.         rigidBod.AddForceAtPosition (fluidDragVecPosX, xpos_face_center);
71.         // LEFT (negX):
72.         Vector3 pointVelNegX = rigidBod.GetPointVelocity (xneg_face_center);
73.         Vector3 fluidDragVecNegX = transform.right.normalized * Vector3.Dot (-transform.right, pointVelNegX) * sa_x * viscosityDrag;
74.         rigidBod.AddForceAtPosition (fluidDragVecNegX, xneg_face_center);
75.  
76.         //--------------DEBUG DEBUG DEBUG:-----------------------------------------------------------------------
77.         // DEBUG VIS positive x,y,z faces centers:
78.         db_xpos.transform.position = xpos_face_center;
79.         db_ypos.transform.position = ypos_face_center;
80.         db_zpos.transform.position = zpos_face_center;
81.  
82.         float debug_multiplier = 1f;
83.         // X Pos debug shape
84.         Vector3 newScalePosX = new Vector3(0.1f, 0.1f, fluidDragVecPosX.magnitude * debug_multiplier);
85.         db_xpos.transform.localScale = newScalePosX;
86.         db_xpos.transform.rotation = Quaternion.LookRotation (fluidDragVecPosX);
87.         // Y Pos debug shape
88.         Vector3 newScalePosY = new Vector3(0.1f, 0.1f, fluidDragVecPosY.magnitude * debug_multiplier);
89.         db_ypos.transform.localScale = newScalePosY;
90.         db_ypos.transform.rotation = Quaternion.LookRotation (fluidDragVecPosY);
91.         // Z Pos debug shape
92.         Vector3 newScalePosZ = new Vector3(0.1f, 0.1f, fluidDragVecPosZ.magnitude * debug_multiplier);
93.         db_zpos.transform.localScale = newScalePosZ;
94.         db_zpos.transform.rotation = Quaternion.LookRotation (fluidDragVecPosZ);
95.  
96.  
97.     }
98. }

-- So I feel like the code is pretty ugly but i'm not sure what kind of Unity functions might be able to work better/faster/more efficiently.

Some Areas i'm unsure of:
-FixedUpdate or different function? Should calling AddForceAtPosition be done during FixedUpdate, or is there a better spot in the main game-loop that would work better?
-This code seems to work ok for linear velocities, but the rigidbody box can spin around its center of mass without any resistance. I'm planning on eventually linking multiple bodies with hingeJoints in the future, which would likely prevent this situation, but should I include an angular drag component, or use the Inertia Tensor?
-I'd like this code to run as fast as possible - are there any obvious things I could do to optimize this?

 

Hi, Eaclou
Interesting topic..

- FixedUpdate is definitely the place to be for AddForceAtPosition. You might want to multiply the force vector in AddForceAtPosition by Time.deltaTime to make sure the code works the same with different fixed time steps.

- It works for linear velocities, but not angular because you are using AddForceAtPosition with a force vector that is always directed exactly towards the centre of mass. That will never give you any angular momentum. Maybe if you used negative pointVelPosX.normalized instead of -transform.right as the force direction?

- About the performance.. transform.right/up/forward are already unit vectors (magnitude = 1), so no need to write transform.right.normalized. That would be basically wasting a Vector3.Distance call each time and that includes calculating a square root.
Also, you might want to cache all those transform.up/right/forwards, I'm pretty sure there is a Quaternion operation going on internally each time you use them (transform.right = transform.rotation * Vector3.right).
Also, if your object is symmetrical, the point velocity at the back side would be the inverse of the point velocity at the front side, so you could do with half as many GetPointVelocity calls.

Cheers,
Pärtel

 

Thanks for the feedback!

I streamlined the code based on your feedback so it should be more lightweight as far as run-time calculations (hard to tell how much since I've only been using <6 shapes at a time, but I know those savings would add up in a more complicated scenario).

I'm going to hold off worrying about angular velocity, since I tested connecting 6 "paddles" together with motor hinge joints (which rotate the blocks from their ends rather than centers of mass) and was able to get a snake-like creature to move pretty similarly to the karl sims video. I can go back and make the physics model more complicated if I have to, but i'm trying to keep it as simple and fast as possible in order to get believable locomotion.

Edit:

Here's a video of my first stab it it, and the code if anyone is interested in doing something similar:

Code (CSharp):

1. using UnityEngine;
2. using System.Collections;
3.  
4. public class sub_box_fluid_drag_script : MonoBehaviour {
5.     public float viscosityDrag = 1f;
6.     private Rigidbody rigidBod;
7.  
8.     // Surface Areas for each pair of faces (neg x will be same as pos x):
9.     private float sa_x;
10.     private float sa_y;
11.     private float sa_z;
12.  
13.     // Use this for initialization
14.     void Start () {
15.         rigidBod = GetComponent<Rigidbody> ();
16.  
17.         // Calculate surface areas for each face:
18.         sa_x = transform.localScale.y * transform.localScale.z;
19.         sa_y = transform.localScale.x * transform.localScale.z;
20.         sa_z = transform.localScale.x * transform.localScale.y;
21.     }
22.  
23.     void FixedUpdate () {
24.         // Cache positive axis vectors:
25.         Vector3 forward = transform.forward;
26.         Vector3 up = transform.up;
27.         Vector3 right = transform.right;
28.         // Find centers of each of box's faces
29.         Vector3 xpos_face_center = (right * transform.localScale.x / 2) + transform.position;
30.         Vector3 ypos_face_center = (up * transform.localScale.y / 2) + transform.position;
31.         Vector3 zpos_face_center = (forward * transform.localScale.z / 2) + transform.position;
32.         Vector3 xneg_face_center = -(right * transform.localScale.x / 2) + transform.position;
33.         Vector3 yneg_face_center = -(up * transform.localScale.y / 2) + transform.position;
34.         Vector3 zneg_face_center = -(forward * transform.localScale.z / 2) + transform.position;
35.  
36.         //=== FOR EACH FACE of rigidbody box: ----------------------------------------
37.         //=== Get Velocity: ---------------------------------------------
38.         //=== Apply Opposing Force ----------------------------------------
39.  
40.         // FRONT (posZ):
41.         Vector3 pointVelPosZ = rigidBod.GetPointVelocity (zpos_face_center); // Get velocity of face's center (doesn't catch torque around center of mass)
42.         Vector3 fluidDragVecPosZ = -forward *     // in the direction opposite the face's normal
43.                                     Vector3.Dot (forward, pointVelPosZ)   // get the proportion of the velocity vector in the direction of face's normal
44.                                     * sa_z * viscosityDrag;   // multiplied by face's surface area, and user-defined multiplier
45.         rigidBod.AddForceAtPosition (fluidDragVecPosZ*2, zpos_face_center);  // Apply force at face's center, in the direction opposite the face normal
46.         // the multiplied by 2 is for the opposite symmetrical face to reduce # of computations
47.  
48.         // TOP (posY):
49.         Vector3 pointVelPosY = rigidBod.GetPointVelocity (ypos_face_center);
50.         Vector3 fluidDragVecPosY = -up * Vector3.Dot (up, pointVelPosY) * sa_y * viscosityDrag;
51.         rigidBod.AddForceAtPosition (fluidDragVecPosY*2, ypos_face_center);
52.        
53.         // RIGHT (posX):
54.         Vector3 pointVelPosX = rigidBod.GetPointVelocity (xpos_face_center);
55.         Vector3 fluidDragVecPosX = -right * Vector3.Dot (right, pointVelPosX) * sa_x * viscosityDrag;
56.         rigidBod.AddForceAtPosition (fluidDragVecPosX*2, xpos_face_center);
57.        
58.     }
59. }
60.  

 

perhaps you should consider developing and refining this into a saleable product.

 

That looks good, eaclou. You're right, the snake moves very much like Karl Sims' evolved snake.

The physics part of your work sounds similar to what I'm doing for one of the forces in my speedboat simulator: http://www.performancesimulations.com/wp/speedboat-simulator/

I do the forces on the individual triangles of the boat mesh. I precompute the center positions and poly areas and store them in local coords, then just transform them into world space from there. Since you're dealing with quads essentially, it makes sense that you're using the quad center instead (half as many forces as two triangles). It looks like you're computing the face positions every step though directly into world coordinates. Maybe you could precompute those locally and then transform them into world space after the fact? That'd be a lot cheaper computationally.

If you're really looking for serious speed, most of your code there can be optimized to be probably 1000% faster or more by getting rid of the math overloads and all the function calls like Dot() and GetPointVelocity() and so forth. Just expand out the computations manually line by line and you might be amazed at the difference.

I wrote an article on it here: http://www.performancesimulations.c...ses-in-unitys-physics-or-any-math-heavy-code/

Caching the transforms like Partel said is a good idea too. If you were to do that, here's an example optimization you could make (not tested):

Code (csharp):

1.  
2. //your original code
3. Vector3 fluidDragVecNegZ = transform.forward.normalized * Vector3.Dot (-transform.forward, pointVelNegZ) * sa_z * viscosityDrag;

Optimized (untested, this is just to give an idea of what I'd do):

Code (csharp):

1.  
2. cachedTransformForward = transform.forward; //Lose the "normalized" like Partel said, it's already normalized)
3.  
4. Vector3 fluidDragVecNegZ;
5.  
6. //Do the dot product first, but break it up and do it manually like this instead of calling the Dot() method:
7. float dot;
8. dot = -cachedTransformForward.x * pointVelNegZ.x +
9.     -cachedTransformForward.y * pointVelNegZ.y +
10.     -cachedTransformForward.z * pointVelNegZ.z;
11.  
12. fluidDragVecNegZ.x = cachedTransformForward.x * dot * sa_z * viscosityDrag;
13. fluidDragVecNegZ.y = cachedTransformForward.y * dot * sa_z * viscosityDrag;
14. fluidDragVecNegZ.z = cachedTransformForward.z * dot * sa_z * viscosityDrag;
15.  

Of course this means that function will get a lot bigger (this is only expanding one of your lines), but it'll run a lot faster.

Depending on how many forces you'll have on each rigid body, I might get rid of the AddForceAtPosition calls too by computing the forces and torques separately on my own. While they're being computed you can add up the total resultant force and torque and add it in a single AddForce and AddTorque call. This is what I do in my boatsim. No matter how many forces I have (probably 2000-3000 per boat right now) there is only one AddForce and one AddTorque call at the end.

If you did all this on the triangles separately, you'd have twice as many forces but they would no longer be centered exactly in the middle of the face. Each quad is two triangles with the center (force application point) offset to the side somewhat. As the quad spins, each triangle center will have a different velocity and force which will add in the damping you're looking for, at least to an extent. It won't be totally physically correct but will probably be good enough for what you're doing.

I'd consider redoing this to use the triangles in the mesh like I do instead if you think this may ever be important. It's a lot more work, but it might be worth it. You might find that it simplifies things when it comes to the actual physics being on a per triangle basis instead of six forces/faces per box. You can get rid of a lot of your x's and y's and z's that way.

It may become easier to add more forces later if you want too. I've got Reynolds number driven skin friction and all kinds of goodies going on in my boat sim on each triangle separately. I also subdivide the triangles at the water surface so a triangle can be partly submerged and partly exposed: http://www.performancesimulations.com/wp/c-waterline-mesh-splitting/

If you're doing virtual creatures, maybe they could evolve to walk on land or turn into birds this way after a few million CPU hours.

 

wow, this was nice reading. i've learned a lot, and bookmarked it for later

 

Thanks for the tips and the information on your site, Todd Wasson.

I never would have realized that the built-in math functions would be slower than directly implementing the basic math functions, but it makes sense in retrospect (I'm coming from the art side of game dev, trying to skill up at coding). Your speedboat sim looks really cool -- I probably won't need to get quite that involved with the physics model, but it still uses some tricks that I think will be helpful when I get further into this project.

What exactly do you mean by:

I was using the transform.forward etc. vectors because it was the only effective way I had found to get the proper vector of the surface normal, and the physics functions (getvelocityatpoint and addForce) seemed to operate in world coordinates so it seemed to work. I'll look into what I could do in local coordinates.

Unfortunately, I'm not going to have much time to work on it for the next month or so due to Work, but i'll post updates when i'm again able to focus some time and energy on it

I'm definitely going to develop it further, but i have a ways to go before it would be polished and presentable enough for me to feel ok charging money for it. I'll be posting on the forums when I can get back into it, but I still have a lot of work to do on the other aspects of it, like the chromosome/crossover model, neural nets, and UI/editors/data visualizations, etc.


Also, does anyone know if there is a simple way to have a Unity program execute its main physics loop "offline" as fast as the processor allows in the background? I can anticipate wanting the ability to churn through the simulations at faster-than-real-time speeds, but am not sure if that would require hand-coding a custom physics system separate from the built-in PhysX?

 

Most people don't seem to realize that and don't seem to check. The problem isn't that the Unity functions are slow per se, it's just that any function call is slower than no function call. With engines structured a bit differently you don't always have that problem with overloads and function calls like these because the compiler will inline the function calls. In effect doing basically the same thing I did in the article. With Unity it can't inline it probably because the Unity engine is a dll while your code is separate from all that. No big deal really, I just consider it good practice. It helps burn the math into my brain.

You might be right. What I had in mind was doing your face position calculations in local space at startup, then getting the world coords by manually multiplying by the object's rotation matrix. Now that I look at it again, I'm not sure which would be faster if both approaches were fully unrolled. In that case it's just a matter of number of operations. Your computations for each face don't have many of them so your way might be faster. I'm too lazy to check right now though, but it may be something to check if you're looking for even more speed. That's of course of vital importance in an evolution simulation like this, so every little bit you can eek out will help.

Good question. I'm not sure, but I'd try increasing the time passage scale (I forget the variable name) and see if it increases the number of FixedUpdate calls per second or if it just ends up changing the fixedDeltaTime in the simulation. I.e., the idea there being to run the simulation in fast forward. If it works (that's a big "if"), you might be able to tune the fixedDeltaTime to be as large as you can get away with in the early stages of evolution when they're struggling to do anything, then maybe decrease it later when the creatures get better. The bigger it is of course the farther away from real time you could get. Then I'd increase the time passage rate to just barely fill up the CPU and bring the frame rate to its knees so you're at maybe 90% physics or something.

Not sure if that will work, but it may be an option. In VRC years ago I was experimenting with neural networks to drive the AI cars in a way that's probably pretty similar to what you're doing. I'd accelerate time in much the same way, effectively by calling my version of FixedUpdate() more often than you're supposed to. If you can find the variable in Unity and double the rate of time passage and also see a doubling in the number of FixedUpdate() calls without Unity modifying your fixedDeltaTime, you could probably do the same thing I did in my evolution simulations: Run the whole simulation in fast forward enough to load the CPU with nearly 100% physics calculations by getting Unity to call FixedUpdate() as many times per second as you can get away with.

By the way, the experiments worked well. My AI was able to learn to drive the cars much more quickly than I could hand tune the controller variables. Evolution is powerful stuff and it's fascinating to watch the systems unfold!

 

By the way, when you get back to it and ever reach a point where you'd like some help optimizing those areas, shoot me a PM with some code. I'd be happy to help.

Another tip: I found particle swarm algorithms for optimizing neural network weights to work really well and be easier to deal with than genetic algos. In your case it might be a departure though because then you wouldn't be modeling the genes quite as literally which might defeat the point of the simulation. Food for thought anyway.