Add combined flow rate boundary condition

What is the requested feature?
In GEOS, we had a combined flow rate boundary condition, where the user would request a total flow be allocated to a set of faces. The solver then would approximate the flow contributions to the faces in the set based on their current permeability. With these flow estimates, the BC was applied using the same approach as the typical flux BC.

This is useful because the user can approximate the behavior of wells without dealing with the additional requirements required by them. I've looked at the old approach, and put together this psuedo-code to show the approach:

# Calculate the average location of faces, total flow rate
n_inlets = 0
sumK_KP = [0.0, 0.0]
qTotal = 0.0
reference_position = [0, 0, 0]
for face in set:
    if not face.ghost:
        reference_position += face.position
        n_inlets += 1
        q_total += bc.get_flow(face, time)
reference_position /= n_inlets

# Estimate flow contributions
for face in set:
    if not face.ghost:
        sumK_KP[0] += face.permeability
        if apply_gravity:
            dx = face.position - reference_position
            gravity_magnitude = dot(dx, gravity_vector)
            sumK_KP[1] += face.permeability * (face.pressure - face.density * gravity_magnitude)
        else:
            sumK_KP[1] += face.permeability * face.pressure

# Choose injection pressure, trial flow rates
q_faces = []
q = 0
if ((n_inlets > 0) && (sumK_KP[0] > 0)):
    sumK_KP[1] += q_total
    injection_pressure = sumK_KP[1] / sumK_KP[0]
    for face in set:
        if use_gravity:
            dx = face.position - reference_position
            gravity_magnitude = dot(dx, gravity_vector)
            q = face.permeability * (injection_pressure - face.pressure + face.density * gravity_magnitude)    
        else:
            q = face.permeability * (injection_pressure - face.pressure)
        q_faces.append(q)
            
# Apply q_faces as the normal mass-rate boundary