Optional valve-state freezing for piecewise valves to prevent Newton cycling with 3D-0D coupling

[dseyler]

Problem

The current piecewise valve implementation allows the valve state to flip between 3D newton iterations, leading to Newton cycling and divergence. For GenBC, we fixed this by freezing the valve state at the beginning of each 3D timestep according to the previous timestep's pressure state. While this does lead to occasional pressure/flow spikes at transitions, this is avoidable with a sufficiently small time step and seems to be the simplest solution to the valve cycling problem.

Solution

I have been implementing the feature in this branch. Let me know your thoughts if there would be a better way to solve this problem.

Additional context

No response

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[aabrown100-git]

Can you describe a little bit how this is implemented? Also how the option is exposed to the user?

[dseyler]

The user sets "freeze_piecewise_valve_state" as a boolean in the simulation_parameters block of the input JSON. It defaults to false if not provided. There's probably a better way to do this, so open to ideas. Maybe exposing it to the user in the valve block would make more sense.

My current implementation in the JSON:

"simulation_parameters": {
  "number_of_cardiac_cycles": 1,
  "number_of_time_pts_per_cardiac_cycle": 689,
  "freeze_piecewise_valve_state": true
}

The parsed flag is then propagated onto the Model object that the valve blocks query via model->freeze_piecewise_valve_state.

When set to true, the valve state is evaluated once per timestep from the previous converged state and the resistance is held constant across all Newton iterations.

A per-step hook is added to Block (src/model/Block.h), defaulting to a no-op. It's fired once per timestep, before the Newton loop, and receives the previous step's converged solution:

virtual void prepare_step(const ...& y_old, const ...& ydot_old);  // default: no-op

Model::prepare_step dispatches the hook to every block (src/model/Model.cpp):

void Model::prepare_step(const ...& y_old, const ...& ydot_old) {
  for (auto block : blocks) {
    block->prepare_step(y_old, ydot_old);
  }
}

PiecewiseValve overrides the hook. It evaluates the valve state against the previous converged pressures and caches the result in R_cached (src/model/PiecewiseValve.cpp):

void PiecewiseValve::prepare_step(const ...& y_old, const ...& /*ydot_old*/) {
  if (!model->freeze_piecewise_valve_state) return;   // feature off -> no-op
  double p_in_old  = y_old[global_var_ids[0]];
  double p_out_old = y_old[global_var_ids[2]];
  R_cached = (p_out_old < p_in_old) ? Rmin : Rmax;
}

update_solution, called once per Newton iteration, reads R_cached if freeze_piecewise_valve_state is set to true, otherwise it uses the per-iteration evaluation(src/model/PiecewiseValve.cpp):

if (model->freeze_piecewise_valve_state) {
  resistance = R_cached;                         // constant across the step
} else {
  resistance = (p_out < p_in) ? Rmin : Rmax;
}

The integration drivers call Model::prepare_step once per timestep, immediately before the integrator step, passing the previous step's state (src/solve/Solver.cpp)

this->model->prepare_step(state.y, state.ydot);  // state = previous timestep
state = integrator.step(state, time);

[mrp089]

Have you tried the smooth tanh valves?

[dseyler]

@mrp089, I tried tanh valves and ran into the following problem: If I made the tanh valves smooth enough to prevent Newton cycling, the valves would "leak" and the chamber pressures would differ from their converged 0D-only values. If I made the valves steep enough to prevent leakage, then they would effectively act like diodes and I would run into the same problems. It's possible that if I tuned the 0D-only model using tanh valves, then that wouldn't have been an issue.

I also tried implementing a line search for each Newton iteration, but this just led to slowly non-descending residuals for the 3D Newton Solver rather than leading to convergence

[mrp089]

Makes sense. We also had some trouble setting up the tanh valve parameters (in a 0D-only model). @KatrinKoesler, can you briefly explain how you tuned the valves in your case (and what final values with units you ended up with)? Happy to compare to the valve freezing!

[aabrown100-git]

I think we should also consider implementing line search in svMultiPhysics to prevent this Newton cycling. According to my understanding of the issues, line search should resolve it. I tested this several years ago in svFSI without success, but I still think it's worth trying in svMultiPhysics.

[KatrinKoesler]

I tuned the valve parameters by running simulations with logarithmically increasing valve R_max (R_min was always set to 0) and chose the smallest R_max value where the resulting pv-loop did not change anymore with a further increase in R_max. I did this first for both valves simultaneously, and then individually for the inlet valve as a smaller R_max value is usually sufficient there (the smaller value improved the convergence behavior a little for me). Then I ran simulations with logarithmically decreasing valve steepness and again chose the largest steepness where the pv-loop would not change anymore (with the steepness being too small convergence was also impaired). I recently used this tuning method on a closed loop system with the spherical ventricle, inlet and outlet tanh valves, and six vessel blocks (with different resistances and capacitances mimicking the systemic circulation) and ended up with the following values:
Inlet valve: "Rmin": 0.0, "Rmax": 1.0e10, "Steepness": 0.1
Outlet valve: "Rmin": 0.0, "Rmax": 1.0e11, "Steepness": 1.0
Unit for the resistances: kg m^-4 s^-1, for the steepness: Pa^-1

[dseyler]

I think we should also consider implementing line search in svMultiPhysics to prevent this Newton cycling. According to my understanding of the issues, line search should resolve it. I tested this several years ago in svFSI without success, but I still think it's worth trying in svMultiPhysics.

I have a branch with it implemented in svMultiphysics and also didn't have much success. It seemed to prevent cycling but didn't necessarily help with convergence.

I think there are two issues at play that lead to divergence in 3D-0D cardiac simulations:

• Valve state flipping: A line search might help with this.
• Isovolumic contraction/relaxation: Line search didn't seem to help. The backtrack line search maxed out at the smallest step size at every non-descending timestep. It's possible that the minimum line search step size was still too big or that the tangent is incorrect and giving a wrong Newton direction. My other idea was that the predictor step is landing too far away from the solution, so the Newton iterations never have a chance to recover. Shrinking the time step seems to help considerably with the isovolumic phase. I wonder how hard it would be to implement time substepping to help with this? I think it could speed up our cardiac simulations a lot since we could increase the timestep during the "easy" parts of the cycle and then adaptively substep during the isovolumic phase.

[dseyler]

I tuned the valve parameters by running simulations with logarithmically increasing valve R_max (R_min was always set to 0) and chose the smallest R_max value where the resulting pv-loop did not change anymore with a further increase in R_max. I did this first for both valves simultaneously, and then individually for the inlet valve as a smaller R_max value is usually sufficient there (the smaller value improved the convergence behavior a little for me). Then I ran simulations with logarithmically decreasing valve steepness and again chose the largest steepness where the pv-loop would not change anymore (with the steepness being too small convergence was also impaired). I recently used this tuning method on a closed loop system with the spherical ventricle, inlet and outlet tanh valves, and six vessel blocks (with different resistances and capacitances mimicking the systemic circulation) and ended up with the following values: Inlet valve: "Rmin": 0.0, "Rmax": 1.0e10, "Steepness": 0.1 Outlet valve: "Rmin": 0.0, "Rmax": 1.0e11, "Steepness": 1.0 Unit for the resistances: kg m^-4 s^-1, for the steepness: Pa^-1

That makes sense. Maybe I could try something more methodical like this for the 3D-0D simulations. I just picked a few different steepness values by hand and couldn't find one that worked so pivoted to just freezing the valves across time steps. I've also found that a smaller Rmax helps with convergence (using Rmax = ~6e9 for both inlet and outlet)