Coupling to 0D solver only supported for first equation

[michelebucelli]

Description

In a number of places, the infrastructure for coupling svMultiPhysics to svZeroDSolver assumes that the coupling is happening for the first equation. This means that e.g. in electromechanics simulations (introduced in #578) the structural mechanics equation needs to be the first. If electrophysiology is the first equation, the 0D solver doesn't run, and no error message or warning is printed to make the user aware of this.

One instance of this is

svMultiPhysics/Code/Source/solver/initialize.cpp

Lines 416 to 422 in 042a155

	for (auto& bc : com_mod.eq[0].bc) {
	// Check for coupled faces (Dir, Neu via cplBC) or Coupled BCs
	if (bc.cplBCptr != -1 || utils::btest(bc.bType, static_cast<int>(consts::BoundaryConditionType::bType_Coupled))) {
	com_mod.cplBC.coupled = true;
	break;
	}
	}

The flag com_mod.cplBC.coupled is set to true only if com_mod.eq[0].bc contains some coupled condition. If the coupled condition is in some other equation, then com_mod.cplBC.coupled remains false, and the 0D solver is silently never called.

A similar behavior is happening in a few other places (I could find similar things in Integrator and set_bc.cpp, but there might be a few more).

Reproduction

From the branch of #578, try running an electromechanics simulation coupled to svZeroDSolver for the structure boundary conditions. The first equation must be CEP, the second struct.

Expected behavior

The coupling with the 0D solver should work independently of the ordering of the equations. The ordering of the equations should affect the order in which they are solved in partitioned schemes, but not the BCs that are applied to them.

Should that be impossible (but I don't see why it should), at least an error message should be raised to notify the user that the 0D solver is not being called despite being specified in the XML file.

Additional context

No response

Code of Conduct

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

[ktbolt]

@michelebucelli There does seem to be an assumed equation ordering for some simulation types. For example. I've seen that if multiple equations are solved and fluid/FSI is involved, the code expects fluid/FSI to be the first equation.

[michelebucelli]

@ktbolt Thanks for the feedback!

I don't know about FSI or other cases, but it seems to me that, at least for the coupling to the 0D solver, there should in principle be no issue to make it work independently of equation ordering. I am trying to work on this downstream of #578. What are your thoughts @javijv4?

For other cases where an ordering is assumed, it might be enough to add an exception for when that ordering is violated, so that the user is at least made aware of this, but I'm not sure what those situations might be (and I haven't encountered any other yet).

[alisonmarsden]

I would like to get @aabrown100-git 's opinion on this - he has the most experience in the group with these BC coupling schemes.

[javijv4]

@michelebucelli, I think, ideally, the user should be able to list the equations they want to solve without considering any ordering restrictions, since these arise from the solver implementation and not from a physics constraint.

That being said, I think we need to at least warn the user if the ordering is not what is expected when doing multiphysics simulations.

Also, maybe this is not easy to implement / a bad idea, but could we do a pre-read of all the equations within the input file, just to check the order and the physics type, so then we can read them in the order we need to within the solver?

Last thing, the way the coupled BCs are implemented and appear in the different places within the code (initialize.cpp, read_files.cpp, set_bc.cpp, etc) is not great, and moving to a more object-oriented implementation could make implementing an order-independent equation reading way easier and clearer. But this is probably another issue.

[aabrown100-git]

I agree, in principle there should be no issue with equation ordering. In svMultiPhysics, the 0D contributions to the linear system are condensed out and the result is simply that residual and tangent terms are added to the coupled equation (either fluid or structure). So 0D coupling is implemented as a solution-dependent BC and shouldn't cause equation ordering issues at all. I think it's just a matter of keeping track which equation is 0D-coupled. We may have to assume that only one equation is coupled. I think that's a fair assumption, and doing so would limit the required code changes.

[ktbolt]

All of these issues will be addressed with the refactoring implementation.

[michelebucelli]

Thanks for your feedback @aabrown100-git @javijv4 @ktbolt!

@ktbolt I agree that the issue will be more cleanly addressed as part of the refactoring, but I also think that a fix within the current framework should be possible and relatively uninvasive, and it will be useful towards electromechanics coupling #523.

I've pushed code addressing this issue to this branch.

The main changes are:

• added a member unsigned int equationIndex to cplBCType, to store the index of the equation that is coupled to the 0D solver;
• modified CoupledBoundaryCondition::compute_flowrates and CoupledBoundaryCondition::compute_pressures so that the equation index is used as an offset when computing integrals; some similar changes are in set_bc.cpp;
• Integrator::step advances the 0D model when cEq == com_mod.cplBC.equationIndex, rather than cEq == 0;
• in initialize, when setting the flag com_mod.cplBc.coupled, we loop through equations instead of just checking equation 0; if different equations specify coupled BCs, an exception is raised, as it seems to me this is currently unsupported.

I still have problems with running a simple 3D-0D electromechanics simulation, but I think that might be due to other issues (e.g. solver tolerance configuration). Once I figure these out I'll open a PR.

Also, maybe this is not easy to implement / a bad idea, but could we do a pre-read of all the equations within the input file, just to check the order and the physics type, so then we can read them in the order we need to within the solver?

@javijv4 I'm not sure what you mean by this: do you mean to enforce an equation order that is independent of what is specified in the XML file (e.g. always solve CEP before struct, regardless of how they appear in the file)?

Last thing, the way the coupled BCs are implemented and appear in the different places within the code (initialize.cpp, read_files.cpp, set_bc.cpp, etc) is not great, and moving to a more object-oriented implementation could make implementing an order-independent equation reading way easier and clearer. But this is probably another issue.

I agree this is a larger issue, but I'd be happy to discuss this later on. In this respect, I noticed there's two classes to store data about coupled conditions, cplBCType and CoupledBoundaryCondition, and I'm not clear on what their intended purposes are and how they differ. Do you have any insight on this?

[javijv4]

@michelebucelli

@javijv4 I'm not sure what you mean by this: do you mean to enforce an equation order that is independent of what is specified in the XML file (e.g. always solve CEP before struct, regardless of how they appear in the file)?

Yes, that is what I meant. The input file will just have a list of physics to be solved, and how we solve them (and in what order) will be totally handled inside the #solver.

I agree this is a larger issue, but I'd be happy to discuss this later on. In this respect, I noticed there's two classes to store data about coupled conditions, cplBCType and CoupledBoundaryCondition, and I'm not clear on what their intended purposes are and how they differ. Do you have any insight on this?

Yes, cplBCType was the legacy implementation. I recently implemented CoupledBoundaryCondition to have a more object-oriented way of dealing with the ZeroD coupling (especially for struct simulations that need a capping surface). I did not want to touch other coupling cases at the time (like RCR), so I left the old cplBCType.
Now that the 1D coupling is implemented, it's probably a good time to give it a go at refactoring all the coupled boundary conditions (ZeroD, OneD, RCR) such that all of them use CoupledBoundaryCondition. I can open an issue if that sounds good.

[ktbolt]

@michelebucelli Thanks for looking into this !

It is good to fix this within the current framework if this can be done cleanly; will allow correct coupling and make it easier to refactor.

[michelebucelli]

Updates on this issue, from this branch. I took the opportunity to address some other small issues related to 3D-0D coupling in the electromechanics setting.

1. The 3D-0D coupling now works correctly independently of equation ordering (as far as I can tell);
2. I had to update the function gnnb (compute area-scaled normal at boundary quadrature points): since the function used the displacement field to compute the normal in deformed configuration, the index of the displacement field in the solution arrays needed to be forwarded to it;
3. I took the opportunity to simplify the gnnb interface a bit, dropping redundant arguments (e.g. no need to pass the size of an array along with that array) and having it return its result, rather than assign it to an output argument;
4. the same changes propagate to those overloads of integ that compute boundary integrals and use gnnb;
5. in many places, both gnnb and integ were used without actually accounting for displacement (this was the case of those simulations not involving structural displacement, for example). To make this cleaner, I introduced overloads to gnnb and integ that do not take the parameters cfg (identifying the mechanical configuration to use) and displacement_index, falling back on the more general case but hiding this from the caller;
6. similar changes apply to the management of caps;
7. the increment for evaluating 0D solver tangent with finite differences was defined as $\epsilon = \max(10^{-5} \text{rms}(Q), 10^{-7})$, with the two constants $10^{-5}$ and $10^{-7}$ hardcoded in calc_der_cpl_bc. The absolute factor $10^{-7}$ is especially problematic because it is dimensional, and so it may not be consistent with the system of units used in the XML file. Both these values can now configured in the XML section <svZeroDSolver_interface>;
8. the output files generated by svZeroDSolver (svZeroD_data, P_svZeroD, Q_svZeroD) are now saved to the same folder as the 3D simulation output (n-procs by default), rather than in the current directory, so that all simulation outputs are bundled together;
9. as requested by @ktbolt, files defining classes in the ActiveStress and IonicModel hierarchies now use upper camel case naming.

I am now able to run a very basic electromechanics + 0D coupled simulation, with an RCR circuit downstream of an idealized ventricle.

I am still having issues with running a slightly more complicated simulation involving valves, with the circuit from Brown et al. (2024), where I observe a sudden failure of convergence in the non-linear solver, leading to the blow up of the solution and subsequent crash. This seems to be related to valve behavior, and is probably independent of the equation ordering problem addressed in this issue, but I am still investigating.

Related to point 8 above, I noticed that the file svZeroD_data is always opened in append mode, so that possible subsequent runs are appended to the same file. Is this done by design, or would it be ok to have the simulation clear the file before writing the header? @javijv4 @aabrown100-git @ktbolt