Advanced Dynamic Simulations

Student Names: Anbang Hu, Yang Yu

Note: Proposal can be found here. When writing the proposal, we did not present in great details what we intended to do, because we had few clue on what was actually coming in this project. As we proceed, the whole picture starts to become clear. Thus, we present detailed tasks in this writeup.

Overview

We have implemented forward Euler and simplectic Euler in A4. In this project, we implement improved Euler and backward Euler for both wave and heat equation within the framework of Scotty3D. When implementing the backward Euler, we try different methods to solve the equation, e.g, Fixed-point iteration, Newton's Method, Solving Linear System. In addition, we support solving the heat equation with forward Euler. (Implementation can be found in mesh.h, mesh.cpp.)

Enhanced Viewer

In order to visualize the effects from advanced solvers, we enhance the viewer from A4 in the following way:

The modification is mostly done in application.h, application.cpp. The detail tasks on temperature visualization will be discribed at the beginning of Heat Equation section.

Wave Equation

Forward Euler

We saw from A4 that forward Euler is not stable. A slight initial offset at the left buttock of the bunny Initial offset

quickly blows up the mesh: Foward Euler (Wave)

Improved Euler

Fortunately, we can improve forward Euler by taking into account of both derivatives at current timestep and next timestep. The formulation for improved Euler for wave equation can be found here.

We code up improved Euler for wave equation in Mesh::improved_euler_wave() in mesh.cpp:

  1. calculate new velocity based on current offset and velocity for each vertex and store them in temporary variables,
  2. calculate new offset based on current offset, current velocity and new velocity for each vertex and store them in temporary variables,
  3. update velocity and offset for each vertex.

Although the mesh still blows up when running for a sufficiently long time, it has better result when running for the same duration as in Forward Euler section: Improved Euler (Wave)

Backward Euler

Improved Euler does not seem to always give stable results. In this section, we look at a more powerful update rule called backward Euler, which is unconditionally stable.

Backward Euler update rule is

u(k+1) = u(k) + tau * u'(k+1)

Since u(k+1) appears on both side of the equation, we need to solve an algebraic equation for u(k+1).

There are numerous ways to solve backward Euler. We experiment on two of them:

Fixed Point Iteration

We describe fixed point iteration for wave equation here and implement it in Mesh::backward_euler_fixedpoint_wave() in mesh.cpp:

  1. create temporary variables to store current velocity and offset for each vertex,
  2. for vertex v, calculate new velocity based on current velocity, v->offset and v->velocity and store it in v->velocity,
  3. for vertex v, calculate new offset based on current offset, v->velocity,
  4. run step 2 and step 3 for several times (in our case, we choose to run 10 times).

Even with large initial offset, the result from fixed point iteration is rather smooth: Backward Euler Fixed Point (Wave)

Newton's method

We descibe Newton's method for wave equation here and implement it in Mesh::backward_euler_newton_wave() in mesh.cpp: 1. create temporary variables to store velocity for each vertex in each iteration, 2. for vertex v, calculate new velocity and new offset using Newton's Method and store it in v->velocity and v->offset (one iteration of Newton Method), 3. run step 2 for several times to get the final result for v->velocity (in our case, we choose to run 10 times), 4. for vertex v, calculate new offset based on current offset, and add v->velocity * timestep.

Newton's method can produce even smoother result: Backward Euler Newton (Wave)

Note: Though stable, both fixed point iteration and Newton's are iterative methods. They need to run several iteration for each vertex update, which slows down the computation, which in turn slows down the video playback.

Heat Equation

To support heat equation, we 1. add a field temperature at each vertex; 2. initialize temperature to 4200 Kelvin at each vertex creation (we follow the physics of blackbody radiation); 3. modify Vertex::laplacian() to allow computation of discrete laplacian for heat equation, 4. follow this post to add temperature_to_rgb() in mesh.cpp to convert temperature into rgb color; 5. modify Mesh::draw_faces() to allow mesh coloring using temperature_to_rgb() based on temperature at each vertex; 6. modify Application::keyboard_event() in application.cpp to allow users to adjust the temperature of selected vertex; 7. modify Application::draw_action() in application.cpp to display temperature of selected vertex.

We use different method to implement the heat equation (forward Euler, improved Euler and backward Euler). After implementing these algorithms, we test them on dae/sky/bunny.dae. At the beginning, all vertices have the same temperature (4200 Kelvin), which will show light red. First we select and increase the temperature at some vertices (to over 30000 Kelvin). The neighborhood of these vertices will become light blue or blue when the temperature keeps increasing. Then we use video playback to observe heat diffusion in the bunny mesh.

Forward Euler

First we use forward Euler to solve the heat equation. We described our algorithm here. We implement it in Mesh::forward_euler_heat() in mesh.cpp: 1. for vertex v, calculate the laplacian and store it in temporary variables, 2. use the formula in the algorihtm to compute the new temperature for each vertex, 3. for vertex v, update the temperature.

The result is shown below. Initially, we select some vertices and increase their temperature (blue regions): Forward Euler (Heat)0 Let the heat diffusion process begin: Forward Euler (Heat)1 Forward Euler (Heat)2 Forward Euler (Heat)3 Forward Euler (Heat)4

Improved Euler

We also use improved Euler to solve the heat equation. We described our algorithm here and code up improved Euler for heat equation in Mesh::improved_euler_heat() in mesh.cpp: 1. calculate new temperature based on current temperature and laplacian for each vertex and store them in temporary variables, 2. calculate new temperature based on current temperature and results in step 1 and store them in temporary variables, 3. update temperature for each vertex.

And the result is shown below. Improved Euler (Heat)0 Improved Euler (Heat)1 Improved Euler (Heat)2 Improved Euler (Heat)3

Improved Euler (Heat)4

Backward Euler

Results from wave equation suggests that backward Euler works better than other update rules. In this section we explore both iterative methods and matrix method:

Fixed Point Iteration

We describe fixed point iteration for heat equation here. We implement it in Mesh::backward_euler_fixedpoint_heat() in mesh.cpp:

  1. create temporary variables to store current temperature for each vertex,
  2. calculate new temperature for each vertex v based on current temperature and v->temperature, and store it in v->temperature,
  3. run step 2 for several times (10 times in our case).

The result is shown below. Backward Euler Fixed Point (Heat)0 Backward Euler Fixed Point (Heat)1 Backward Euler Fixed Point (Heat)2 Backward Euler Fixed Point (Heat)3 Backward Euler Fixed Point (Heat)4

Newton's Method

We describe Newton's method for heat equation here. We implement it in Mesh::backward_euler_newton_heat() in mesh.cpp: 1. create temporary variables to store current temperature for each vertex, 2. calculate new temperature using Newton Method for each vertex v based on current temperature and v->temperature, and store it in v->temperature, 3. run step 2 for several times to get the final temperature for each vertex (10 times in our case).

The result is shown below. Backward Euler Newton (Heat)0 Backward Euler Newton (Heat)1 Backward Euler Newton (Heat)2 Backward Euler Newton (Heat)3 Backward Euler Newton (Heat)4

Power of Matrix

As mentioned in Wave Equation section, backward Euler implemented using an iterative method tends to be fairly slow. One way to speed it up is to assemble everything into a matrix and solve a big linear system.

We describe backward Euler with matrix method for heat equation here. We follow the idea in the section titled "Meshes and Matrices" in these course notes as well as linear algebra package Eigen to implement it in Mesh::backward_euler_matrix_heat() in mesh.cpp:

  1. Assign unique index to the vertices in the mesh;
  2. Build the Laplace and mass matrice (sparse matrix) according to the indices of the vertices;
  3. Solve the linear system that describes the backward Euler update rule.

The result is shown below. Backward Euler Matrix (Heat)0 Backward Euler Matrix (Heat)1 Backward Euler Matrix (Heat)2 Backward Euler Matrix (Heat)3 Backward Euler Matrix (Heat)4

We expected that solving linear systems would be faster than iterative methods. However, it turns out that the opposite is true. A reason may be that iterative methods run for too few iterations (only 10 times). Another reason is that optimization in solving linear systems can be done to accelerate the computation.

Conclusion

In this project, we implement advanced solvers for both wave equation and heat equation. We see that improved Euler gives a slightly more stable result than forward Euler, whereas backward Euler produces the most smooth result. In theory, iterative approximation for backward Euler is slow; so whenever possible, it is faster to implement backward Euler by first converting the problem into a linear system represented as sparse matrix and then solve it. However, our experiments show that solving linear systems may not be as fast as iterating for a few iterations.

Summary of Wave & Heat Equation Solvers We Implement

forward_euler:
1. forward_euler_wave
2. forward_euler_heat

improved_euler:
1. improved_euler_wave
2. improved_euler_heat

backward_euler:
1. backward_euler_wave
    a. backward_euler_fixedpoint_wave
    b. backward_euler_newton_wave (default)
2. backward_euler_heat
    a. backward_euler_fixedpoint_heat
    b. backward_euler_newton_heat (default)
    c. backward_euler_matrix_heat

The default for backward_euler_wave is to use backward_euler_newton_wave. In order to see the effect of backward_euler_fixedpoint_wave, select the corresponding option in backward_euler_wave and recompile the code. In order to see the effect of backward_euler_fixedpoint_heat or backward_euler_matrix_heat, select the corresponding option in backward_euler_heat and recompile the code.

Key Commands

Command Key
Enter animate mode A
Enter wave action W
Change simulation to use Forward Euler F
Change simulation to use Backward Euler B
Change simulation to use Improved Euler Y
Change simulation to use Symplectic Euler S
Increase thermal diffusivity for Physical Simulation SHIFT + .
Decrease thermal diffusivity for Physical Simulation SHIFT + ,
Increase temperature at selected vertex .
Decrease temperature at selected vertex ,

To play heat diffusion simulation, press SPACE in Wave action. The other key commands are same in the user guide in Scotty3D.

References

We list all the references we have used in the following: