Make SAP the default discrete solver

[amcastro-tri]

This PR effectively makes SAP our default contact solver.
Towards #19322.
This does not address #21225 just yet (though we could).

Though not strictly necessary, I believe we'd want to update all uses of set_discrete_contact_approximation() to set our recommended default. I did this with the deformable bubble gripper sim as a proof of concept.

We mark this as a breaking change, since changing the underlying contact approximation can lead to slightly different simulation results.
This might not be noticeable in most cases, but users should at least be aware of this.
In particular, most users of SAP ommit setting disspation values and use the defult value of relaxation_time, which is non-zero.
However, the newly adopted approximation uses the more physical (and desirable) Hunt & Crossley model, with a default value of dissipation equal to zero.
Therefore users must be aware that from now on, modeling inelastic collisions requires explicitly setting a value for the hunt and crossley dissipation, drake:hunt_crossley_dissipation when setting proximity properties through a model file.

---

This change is 

[sherm1]

@amcastro-tri consider making this a Draft PR until it's ready for review

[amcastro-tri]

Reviewable status: 1 unresolved discussion, needs platform reviewer assigned, needs at least two assigned reviewers, labeled "do not merge"

---

-- commits line 2 at r2:
Note to self:
In the next revision update comment to something that includes why this is a breaking change.

"This PR changes the default discrete solver from TAMSI to SAP, using the kSimilar contact approximation.

The update makes MultibodyPlantConfig::discrete_contact_approximation="similar".

We mark this as a breaking change, since changing the underlying contact approximation can lead to slightly different simulation results.
This might not be noticeable in most cases, but users should at least be aware of this.
In particular, most users of SAP ommit setting disspation values and use the defult value of relaxation_time, which is non-zero.
However, the newly adopted approximation uses the more physical (and desirable) Hunt & Crossley model, with a default value of dissipation equal to zero.
Therefore users must be aware that from now on, modeling inelastic collisions requires explicitly setting a value for the hunt and crossley dissipation, drake:hunt_crossley_dissipation when setting proximity properties through a model file."

[amcastro-tri]

+@jwnimmer-tri for feature review please

Reviewable status: 1 unresolved discussion, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers

[jwnimmer-tri]

Anzu CI is happy with this PR -- just the expected two test failures in the change-detector tests for scenario defaults, which we can easily fix when we bump Drake in Anzu.

Your recent note about zero H&C dissipation give me pause and I wonder if we should fix those defaults prior to (or concurrently with) this PR, for a smoother transition. I will give it some thought.

Reviewed 1 of 3 files at r1, 1 of 2 files at r2, all commit messages.
Reviewable status: 4 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

---

multibody/plant/multibody_plant.cc line 299 at r2 (raw file):

  // By default, MultibodyPlantConfig::discrete_contact_approximation and
  // MultibodyPlantConfig::discrete_contact_solver are empty, indicating that
  // TAMSI is the default approximation and solver.

nit This comment is stale now.

---

multibody/plant/multibody_plant.cc line 304 at r2 (raw file):

  DRAKE_DEMAND(MultibodyPlantConfig{}.discrete_contact_solver == "");
  DRAKE_DEMAND(MultibodyPlantConfig{}.discrete_contact_approximation ==
               "similar");

This says "similar" but two lines prior (L301) says "tamsi". They are supposed to be identical, that's the whole purpose of this cross-check.

The MbP member field default in the header needs to change as does L301.

---

multibody/plant/multibody_plant_config.h line 67 at r2 (raw file):

  /// The underlying solver will be inferred automatically. The deprecated code
  /// will be removed from Drake on or after 2024-04-01.
  std::string discrete_contact_solver{""};

To me, it seems strictly better to remove this old field in this PR. Now that the other field has a default, any non-empty value here could easily cause a crash (b/c both cannot be set at once), so there is very little if any point to still keeping this.

Is there any reason I'm missing to keep this around still? If not, please remove it in this PR and clean up the related comments below and throughout dated 2024-04-01. (Remove MbP::set_discrete_contact_solver, etc.)

[amcastro-tri]

Makes sense. Let me know if you need any help with specific Anzu sims.
I deprecated discrete_conact_solver in this revision, though I might need to revisit some tests. I'll take a look at CI later today or early tomorrow.

Reviewable status: 3 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

---

multibody/plant/multibody_plant.cc line 304 at r2 (raw file):

Previously, jwnimmer-tri (Jeremy Nimmer) wrote…

This says "similar" but two lines prior (L301) says "tamsi". They are supposed to be identical, that's the whole purpose of this cross-check.

The MbP member field default in the header needs to change as does L301.

ah, obviously you are right. Fixed.

---

multibody/plant/multibody_plant_config.h line 67 at r2 (raw file):

Previously, jwnimmer-tri (Jeremy Nimmer) wrote…

To me, it seems strictly better to remove this old field in this PR. Now that the other field has a default, any non-empty value here could easily cause a crash (b/c both cannot be set at once), so there is very little if any point to still keeping this.

Is there any reason I'm missing to keep this around still? If not, please remove it in this PR and clean up the related comments below and throughout dated 2024-04-01. (Remove MbP::set_discrete_contact_solver, etc.)

No reason to delay this, I just wanted to get some quick feedback first.
I elliminated the deprecated code in this revision.

[jwnimmer-tri]

Thanks, I'll check Anzu again tomorrow. The code changes here at least LGTM.

Reviewed 7 of 11 files at r3, all commit messages.
Reviewable status: LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

---

multibody/plant/multibody_plant_config.h line 67 at r2 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

No reason to delay this, I just wanted to get some quick feedback first.
I elliminated the deprecated code in this revision.

Yup, figured. Thanks.

[jwnimmer-tri]

Now that the class default is actually changed in the latest revision, Anzu CI shows the same kinds of errors as Drake -- code that needs AutoDiff crashes now when we select SAP instead of TAMSI.

[amcastro-tri]

Here in Drake I see two types of failure:

1. Tests originally setup for TAMSI, with tolerances calibrated for TAMSI, now fail with SAP. Easy fix, make those specific tests run with TAMSI, or updates tolerances.
2. I see contact_visualizer_test needs scalar conversion. Here again we can make that specific test run with TAMSI.

What about the scalar conversion in Anzu? is that a required functionality or could we use TAMSI for those specific tests?

Reviewable status: LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge"

[jwnimmer-tri]

I see contact_visualizer_test needs scalar conversion. Here again we can make that specific test run with TAMSI.

That's just a band-aid. The goal is that users who just want default values don't need to copy-paste a bunch of TAMSI boilerplate everywhere. I will look at the details more closely next week. My current impression is that there is no MbP config setting that is universally best, so we probably need to fix some MbP bugs / missing features before we can switch to SAP by default.

[amcastro-tri]

Let me summarize the bahavior in case it helps.

• For current users of TAMSI, the same defaults that apply to kTamsi apply to kSimilar. The only caveat here is that AutoDiff is not possible through SAP. I see this as a lesser problem right now, since users who want gradients today, do make use of TAMSI purposedly (though now they'll have to say it explicitly, hence a breaking change).
• I believe users who created models for SAP specifically, with no need to specify Hunt-Crossley dissipation (H&C), will probably be the most affected. Here the physical model of dissipation changes, and thus we'd expect small differences. However the biggest difference has to do with the zero/non-zero default value of the dissipation parameter:
  • For point contact: both models have a default, non-zero, value of dissipation. For SAP, a hard coded relaxation_time of 0.1 s and for H&C, the TAMSI default estimated from the penetration allowance.
  • For hydroelastic: SAP still uses the hard coded relaxation_time of 0.1 s while for H&C we now use a value of zero.

The reason we chose a non-zero default value of relaxation_time is because the dissipation constant for kSap is difficult to tune; low values might lead to unexpected energy gains and high values greatly affect the "gliding during slip". Therefore we decided to provide a "good" value by default for users.

For H&C this is not a problem (and therefore the research effort to convexity this model). We can safely provide zero values (with the much improved numerics of the model) and very high values won't give users headaches.

Reviewable status: LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

[amcastro-tri]

I fixed (C++) tests that fail due to tolerance issues.
For tests like the ones on plant_test.py it'd seem that what is missing is support for T != double at least for when there is no contact. We should be able to support that, though I'd PR separately.
Let me know if you agree with this assessment @jwnimmer-tri

Reviewable status: LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

[amcastro-tri]

Reviewable status: 1 unresolved discussion, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

---

multibody/meshcat/test/contact_visualizer_test.cc line 42 at r5 (raw file):

    // we want to unit test scalar conversion for a diagram containing a contact
    // visualizer, we need a scalar convertible plant. Therefore we use the
    // kTamsi approximation.

note to self: I'd revert this if we can support AutoDiff when there is no contact with SAP

[jwnimmer-tri]

Thanks, I still need to do some homework here before I reply in full. Hopefully I can do that tomorrow, but we'll see what the week brings.

I'm not as worried about the physical constant defaults. Instead, what I want to do is go through all of the MbP options (the 4 contact approximation options plus continuous time) crossed with scalar types (all 3) and take an inventory of which features each one supports correctly, supports incorrectly (prints a warning and computes wrong results), and doesn't support (throws).

[amcastro-tri]

I like that. Notice you can reduce the set by considering solvers (MbP::get_discrete_contact_solver()) instead of approximations (MbP::set_discrete_contact_approximation()). The former only has two outcomes (kTamsi and kSap).

I just went through the code and I am pretty certain this is the set of supported/unsupported entries in that table. Hopefully this helps.

• DiscreteContactSolver::kTamsi: In the abscence of contact, all scalar types are supported. When there is contact, the scalar support is constrainted by the geometry scalar support. We get an exception if we perform a geometry query with a scalar that is not supported.
• DiscreteContactSolver::kSap: Symbolic is not supported for discrete updates, MbP instrospection queries and state depndent calcs (e.g. mass matrix) are supported. AutoDiff is supported only if there are no constraints.

Reviewable status: 1 unresolved discussion, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

[jwnimmer-tri]

I've created the start of my assessment chart here (TRI only). It's still very early, but I'm working on it.

Big picture, as of today neither TAMSI nor SAP is good enough to be a safe default. Since we don't want to invest anymore in TAMSI moving forward (e.g., by adding support for implicit PD to it), that means the path forward for this PR will be to fix what's broken/missing in SAP (in separate PRs) until the feature matrix is good enough, at which point we can land this tiny & trivial PR.

Reviewed 1 of 11 files at r3, all commit messages.
Reviewable status: 5 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

a discussion (no related file):
Working (physical constants)

Once everything else is resolved, we'll need to circle back on the question of default physical constants. It seems like it will be a few careful cross-checks and not a big risk, so I'm deferring poking at it for now.

It's also possible that either #20764 or something like it could (or will have) already merge ahead of this PR, and already have provided sensible default values -- in which all this PR review needs is to cross-check, and probably not change anything further.

---

a discussion (no related file):
Working (AutoDiff)

As best I can tell so far, we absolutely must not merge this PR until we have improved the AutoDiff story for SAP. The fallout seems pretty severe. I'll be working on a the smallest possible "minimum list" of must-have features before we can enable SAP, but suffice to say that any throwing from MbP/SAP that isn't due to e.g. SG throwing is at least very suspicious, and probably will be a blocker.

---

a discussion (no related file):
Working (symbolic)

We can probably entertain some loss of default symbolic expression features (symbolic users can still opt-out back to TAMSI or whatever they use now), but I still want to be sure I have my head around exactly what we're gaining/losing. Ideally we'd have a tabular support matrix in the docs.

---

---

examples/hydroelastic/python_ball_paddle/contact_sim_demo.py line 56 at r5 (raw file):

    plant.Finalize()

    meshcat = Meshcat()

Because Meshcat should nearly always be a singleton, it should nearly always be an input argument, not a fresh return value. The outermost main function should create a Meshcat() and pass it in here.

Anyway, since this appears to be a debugging aid, not a functional change related to the solver selection, please open this up as a separate PR to avoid cluttering this one.

---

multibody/meshcat/test/contact_visualizer_test.cc line 42 at r5 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

note to self: I'd revert this if we can support AutoDiff when there is no contact with SAP

This is a signal that MbP+SAP is not sufficient for users yet. The workflow exercised by this test is a not-uncommon operation, that we must support with our default plant configuration.

[amcastro-tri]

Thank you so much @jwnimmer-tri, that table is super useful. I added a few comments/questions.

I agree with you on the path forward. I'll create a separate PR for each of those and work them out one at a time.

Reviewable status: 5 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge"

[amcastro-tri]

Reviewable status: 5 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge"

a discussion (no related file):

Previously, jwnimmer-tri (Jeremy Nimmer) wrote…

Working (AutoDiff)

As best I can tell so far, we absolutely must not merge this PR until we have improved the AutoDiff story for SAP. The fallout seems pretty severe. I'll be working on a the smallest possible "minimum list" of must-have features before we can enable SAP, but suffice to say that any throwing from MbP/SAP that isn't due to e.g. SG throwing is at least very suspicious, and probably will be a blocker.

Working on #21370 this week. PTAL the issue @jwnimmer-tri, I believe that captures the state of our discussion well. Thanks.

---

[jwnimmer-tri]

Reviewable status: 5 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

a discussion (no related file):

I believe that captures the state of our discussion well.

Yes, looks great, thanks.

---

[amcastro-tri]

Reviewable status: 5 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge"

a discussion (no related file):

Previously, jwnimmer-tri (Jeremy Nimmer) wrote…

Working (symbolic)

We can probably entertain some loss of default symbolic expression features (symbolic users can still opt-out back to TAMSI or whatever they use now), but I still want to be sure I have my head around exactly what we're gaining/losing. Ideally we'd have a tabular support matrix in the docs.

TamsiSolver support of symoblic is a lie. What really happens is that the solver can only handle the case for zero contacts. Instead of compiling TamsiSolver on symobolic, the case for zero contacts should be handled separately in the driver and never instantiante a solver on symoblic::Expresion, it makes no sense.
In practice, users would get an exception from SceneGraph if and only if there is contact (again, for T = symbolic::Expression).

For SAP, we should enable the case for zero contacts in the same way we have it for TAMSI. That is, we should allow to perform discrete updates when there are no constraints. Again, the branch for zero constraints would be dealt with in the driver, not the solver.

Solvers should never be even compiled for symbolic::Expression.

---

[amcastro-tri]

Reviewable status: 5 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge"

a discussion (no related file):

Previously, jwnimmer-tri (Jeremy Nimmer) wrote…

Working (physical constants)

Once everything else is resolved, we'll need to circle back on the question of default physical constants. It seems like it will be a few careful cross-checks and not a big risk, so I'm deferring poking at it for now.

It's also possible that either #20764 or something like it could (or will have) already merge ahead of this PR, and already have provided sensible default values -- in which all this PR review needs is to cross-check, and probably not change anything further.

@rpoyner-tri, where is currently documented what we do (or will do) with physical constants for contact?

When using Hunt&Crossley model with regularized friction (point or hydro), all supported approximations use the same physical parameters (friction, Hunt&Crossley disspation, point stiffnes or hydro modulus, and stiction tolerance).

The only one approximation that is different is kSap, which uses a linear model of dissipation (parameterized by relaxation_time) and an intrinsic regularization of friction (parameterization not exposed. stiction tolerance is ignored).

---

[jwnimmer-tri]

Reviewable status: 5 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

a discussion (no related file):
To be clear: "support T=Expression" does not imply "compute the analytical forms using unbound variables for q,v". A completely valid use case for T=Expression is to have all of the concrete Expression values seen by MbP filled with plain floating point values, in which case the results should look like T=double. (Other systems in the diagram can still have symbolic forms using variables, even if MbP doesn't.)

It's fine if we have a support matrix with certain things marked as "not supported", but I want to make sure we're speaking the same language here. The question is "what can/will we do in the code when T=Expression" which is not the same as "compute symbolic dynamics".

If we call ExtractDoubleOrThrow and it doesn't throw, then we can use the same logic thereafter as for T=double.

Solvers should never be even compiled for symbolic::Expression.

Nothing in my OP suggested that they should be. It's fine by me to short-circuit / re-route prior to reaching the solver.

---

[jwnimmer-tri]

Reviewable status: 5 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

a discussion (no related file):

Previously, amcastro-tri (Alejandro Castro) wrote…

@rpoyner-tri, where is currently documented what we do (or will do) with physical constants for contact?

When using Hunt&Crossley model with regularized friction (point or hydro), all supported approximations use the same physical parameters (friction, Hunt&Crossley disspation, point stiffnes or hydro modulus, and stiction tolerance).

The only one approximation that is different is kSap, which uses a linear model of dissipation (parameterized by relaxation_time) and an intrinsic regularization of friction (parameterization not exposed. stiction tolerance is ignored).

The pull request for default proximity property config is #21366, in case that's what you mean?

---

[amcastro-tri]

Reviewable status: 5 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge"

a discussion (no related file):

Previously, jwnimmer-tri (Jeremy Nimmer) wrote…

The pull request for default proximity property config is #21366, in case that's what you mean?

Thanks, I just didn't know where to start looking.
What I need is in DefaultProximityProperties (geometry/scene_graph_config.h).
Thus far I see that both relaxation_time and hunt_crossley_dissipation are both nullopt, which I believe today it means we keep using hard-coded defaults.

Today the default values are relaxation_time=0.1 s for both point and hydro.
For H&C dissipation, it is non-zero (estimated by MbP based off penetration_allowance) and zero for hydroelastics.

IMO, it'd be nice they both had non-zero default values, regardless of the model used. E.g. relaxation_time=0.1 s and hunt_crossley_dissipation = 20 s/m.

---

[amcastro-tri]

Reviewable status: 5 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge"

a discussion (no related file):

Previously, jwnimmer-tri (Jeremy Nimmer) wrote…

I believe that captures the state of our discussion well.

Yes, looks great, thanks.

#21431 Allows to use AutoDiffXd with SAP.
#21441 confirms this.

---

a discussion (no related file):

Previously, jwnimmer-tri (Jeremy Nimmer) wrote…

To be clear: "support T=Expression" does not imply "compute the analytical forms using unbound variables for q,v". A completely valid use case for T=Expression is to have all of the concrete Expression values seen by MbP filled with plain floating point values, in which case the results should look like T=double. (Other systems in the diagram can still have symbolic forms using variables, even if MbP doesn't.)

It's fine if we have a support matrix with certain things marked as "not supported", but I want to make sure we're speaking the same language here. The question is "what can/will we do in the code when T=Expression" which is not the same as "compute symbolic dynamics".

If we call ExtractDoubleOrThrow and it doesn't throw, then we can use the same logic thereafter as for T=double.

Solvers should never be even compiled for symbolic::Expression.

Nothing in my OP suggested that they should be. It's fine by me to short-circuit / re-route prior to reaching the solver.

#21441 Reveals that wether there is geometry or not, or wether the multibody system is in a state of contact or not, both TAMSI and SAP throw when attempting to perform a discrete update.

---

[jwnimmer-tri]

Reviewable status: 5 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers, labeled "do not merge" (waiting on @amcastro-tri)

a discussion (no related file):

Previously, amcastro-tri (Alejandro Castro) wrote…

#21431 Allows to use AutoDiffXd with SAP.
#21441 confirms this.

Reminder: there are some code changes below (that were in response to the prior lack of autodiff) that need to be reverted now.

---