I deal with a similar situation where vacuum valves prevent negative pressures. Using the toolkit, I perform an initial solution, then select the node with highest hgl where pressure is less than 0 (i.e. node closest to source with psi < 0). I call a function to artificially raise the terminating reservoir elevation until the pressure at my selected node = 0. This gives me reasonably good results between the source reservoir and selected node. I export these results. For the remainder of the system I disable the source reservoir and follow the same procedure, assuming a 0 psi reservoir at my selected node. Continue until the solution has all pressures >=0. It isn’t perfect, but provides results close enough for what I need.