{
  // force residual is the net force on the model
  // this should got to zero in a converged steady state model
  // should be altered for parallel runs
  vector netForce = vector::zero;
  forAll(mesh.boundary(), patchi) {
    netForce +=
    sum(
        mesh.Sf().boundaryField()[patchi]
        &
        (
            2*mu.boundaryField()[patchi]*symm(gradU.boundaryField()[patchi])
            + lambda*tr(gradU.boundaryField()[patchi])*I
        )
    );
  }
  forceResidual = mag(netForce);
}
