I'm considering implementing (just for simplicity) the unconstrained implicit optimization based integration for Material Point Method as described in Chenfanfu Jiang's thesis on MPM (the minimization algorithm starts on page 101 – "A.3.1 Unconstrained minimization").

I have difficulties understanding how would be the algorithm implemented.

Are the statements below (in)correct?

The goal of the algorithm is obtaining new grid velocities in response to external and internal forces of the material represented by material points.

The function to be minimized is $E(\mathbf{v}) = \sum_i\frac{1}{2}m_{\mathbf{i}} \lVert \mathbf v_i - \mathbf {v_i}^n \rVert + \Phi(\mathbf {x_i}^n + \Delta t\mathbf v_i)$. (from Jiang's thesis)

It's derivative is $\mathbf h(\mathbf v) = \mathbf M\mathbf v - \Delta t \mathbf f(\mathbf x^n + \Delta t \mathbf v) - \mathbf M\mathbf v^n$.

Based on the thesis I assume the approach to minimizing $E$ is by the Newton method. The Newton method requires inverse of hessian $\mathbf{\mathit{H}}$. However we're in fact interested in $\Delta \mathbf v$ that would minimize (the linearization of) $E$, so we can instead solve $\mathbf{\mathit{H}}\Delta \mathbf v=\mathbf h$.

When solving the system via Conjugate gradient, we only need to know the product $\mathbf{\mathit{H}}\Delta \mathbf v$, which can be thought of as differential $\delta \mathbf h = \mathbf M\Delta \mathbf v - \Delta t^2 \delta \mathbf f(\mathbf x^n + \Delta t \mathbf v) = \mathbf{\mathit{H}}\Delta \mathbf v$.

From [Stomakhin et al. MPM snow paper] we know how to compute $\delta\mathbf{f_i}$ (page 6 – "6 Stress-based forces and linearization") – i.e. force differential for the grid velocity at node $\mathbf i$.

My assumptions:

  1. When solving $\mathbf{\mathit{H}}\Delta \mathbf v=\mathbf h$ via Conjugate Gradient (step 4. in the algorithm)
    1. ${\Delta \mathbf v}$ is a column vector consisting of component-wise grid velocity differentials – $[{\Delta v}_{1x}, {\Delta v}_{1y}, {\Delta v}_{1z},\,\,\, {\Delta v}_{2x}, {\Delta v}_{2y}, {\Delta v}_{2z}, \,\,\,\ldots]$ (length of $\Delta \mathbf v$ is therefore $3n$, where $n$ is the number of grid nodes (containing velocity))
    2. When multiplying $\mathbf{\mathit{H}}\Delta \mathbf v$, we're in fact always processing 3 rows at the same time - i.e. we process one grid node at time (which consists of x,y,z velocity components).
    3. When computing residual $\mathbf r = \mathbf h - \mathbf{\mathit{H}}\Delta \mathbf v$, we can compute its consecutive 3 components corresponding to grid node $\mathbf i$ (like in assumption 1.1.) by writing $\mathbf{r_i} = \left(m_{\mathbf i}\mathbf v_{\mathbf i} - \Delta t \mathbf f_{\mathbf i}(\mathbf x^n_{\mathbf i} + \Delta t \mathbf v_{\mathbf i}) - m_{\mathbf i}\mathbf {v_{\mathbf i}}^n\right) - \left(m_{\mathbf i}\Delta\mathbf v_{\mathbf i} - \Delta t^2 \delta \mathbf {f_i}(\mathbf x^n_{\mathbf i} + \Delta t \mathbf v_{\mathbf i})\right)$.
    4. When computing $\mathbf r \cdot \mathbf r$, I ignore the fact that each 3 consecutive elements of vector $\mathbf r$ correspond to 1 grid node, and just compute the dot product like with any other vector.

Additional question:

  1. Is the line search (step 7 of the algorithm) necessary?
  • 1
    $\begingroup$ Line searches are never necessary unless you want the algorithm to converge to a solution. $\endgroup$ Jan 28 '17 at 13:18
  • $\begingroup$ Thank you @MarkL.Stone. Does this mean that without line search I would not be able to get anywhere near the solution - i.e. line search is not just for speeding the optimization up, but it's needed in order to find meaningful approximations of solutions? $\endgroup$
    – sarasvati
    Jan 28 '17 at 13:31
  • 1
    $\begingroup$ Without line search, the algorithm may or may nor converge. So use line search. It's possible that without line search, you can be lucky and the algorithm converges to a solution. $\endgroup$ Jan 28 '17 at 14:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.