Given a quantum circuit $C_1$ that generates a state $\vert\psi\rangle$ and another circuit $C_2$ that generates $\vert\phi\rangle$, is there a way to construct a circuit that outputs
$$\frac{1}{\sqrt{2}}(\vert \psi\rangle +\vert\phi\rangle)$$
using $C_1$ and $C_2$ as black boxes?
Here's one way to do it. Let's start with some assumptions: here, I'm assuming your circuits $C_1$ and $C_2$ use the same number of qubits. In the drawing, I have used four qubits to illustrate the concept, but that doesn't matter. The answer below does not care about the number of qubits (which I call $n$), just that the two circuits have the same number of qubits.
I also assume that since we have these two circuits $C_1$ and $C_2$, we can create controlled versions of them. This is not hard to do. [edit after @DaftWullie's comment: the question above starts, "Given a quantum circuit..." so I am assuming here that I know what the circuits are. In the general case of completely unknown circuits, this may not be as easy.]
In addition to the number of qubits needed for your circuits, I also need one more qubit to do this work. It's called ctrl$_0$ in the figure.
Start with all the qubits reset to the state $\vert0\rangle$. Then, apply a Hadamard gate to ctrl$_0$, and you have the combined state
$$ \frac{1}{\sqrt{2}}\left(\vert0\rangle + \vert1\rangle\right)\otimes \vert0\rangle_n $$ where the ordering is apparent from above.
Next, apply the controlled-$U_1$ box, which really just contains the circuit $C_1$ that you have given me. The combined state is now
$$ \frac{1}{\sqrt{2}}\left(\vert0\rangle\vert0\rangle_n + \vert1\rangle\right\vert\psi_1\rangle_n) $$
After applying the $X$ gate, we have
$$ \frac{1}{\sqrt{2}}\left(\vert1\rangle\vert0\rangle_n + \vert0\rangle\right\vert\psi_1\rangle_n) $$
And so, after applying the controlled-$U_2$ box, which really just contains the circuit $C_2$ that you have given me
$$ \frac{1}{\sqrt{2}}\left(\vert1\rangle\vert\psi_2\rangle_n + \vert0\rangle\right\vert\psi_1\rangle_n) $$
After applying the Hadamard gate again on ctrl$_0$, this becomes
$$ \left[\frac{1}{\sqrt{2}}\vert0\rangle\otimes\frac{1}{\sqrt{2}}\left(\vert\psi_1\rangle_n + \vert\psi_2\rangle_n\right)\right] + \left[\frac{1}{\sqrt{2}}\vert1\rangle\otimes\frac{1}{\sqrt{2}}\left(\vert\psi_1\rangle_n - \vert\psi_2\rangle_n\right)\right] $$
Now, if you measured the qubit ctrl$_0$, you will either get 0 or 1. If you measured 0, then you know that the $n$ qubits are now in the state
$$\frac{1}{\sqrt{2}}\left(\vert\psi_1\rangle_n + \vert\psi_2\rangle_n\right)$$
If you measured 1, then you know that the $n$ qubits are now in the state
$$\frac{1}{\sqrt{2}}\left(\vert\psi_1\rangle_n - \vert\psi_2\rangle_n\right)$$
Each of these two outcomes has a probability of 1/2. That means if you really only want the state
$$\frac{1}{\sqrt{2}}\left(\vert\psi_1\rangle_n + \vert\psi_2\rangle_n\right)$$
then this circuit will give it to you half the time. The other half the time, you know (because you measured 1) that you need to restart the circuit in hopes of getting the outcome that you'd like. Since the probability is 1/2, sooner or later this ought to happen.
This method is probabilistic, of course, but has the nice benefit of giving you a measurement outcome (0 or 1) that tells you if you've done what you wanted or not.
As Davit Khachatryan's answer points out, the task is impossible / ill-defined, since the desired target state is generally not normalized and it depends on the relative global phases of the two initial states. However, it is possible to rephrase the question so that is meaningful and has an interesting answer.
The two problems -- sensitivity to the global phase and normalization -- can be simultaneously dealt with by rephrasing the problem in terms of density matrices instead of pure states. Indeed, the density matrix $|\psi\rangle\langle\psi|$ is not sensitive to the arbitrary global phase of $|\psi\rangle$.
To make sense of addition, note that the desired target state $\alpha |\psi\rangle + \beta |\phi\rangle$, for any choice of $\alpha, \beta \in \mathbb C$, lies in the 2-dimensional subspace spanned by $|\psi\rangle$ and $|\phi\rangle$. Hence, without loss of generality, we can restrict the problem to this subspace and think of $|\psi\rangle$ and $|\phi\rangle$ as qubit states regardless of their original dimension.
Now, consider a more general problem where you want to somehow continuously interpolate between $|\psi\rangle$ and $|\phi\rangle$. That is, you want the desired target state to lie on a path (on the Bloch sphere) connecting the two original states. Say, you could choose something like $\sin(t) |\psi\rangle + \cos(t) |\phi\rangle$, appropriately normalized. While this seems like a reasonable choice, the path can vary wildly depending on the global phase of $|\psi\rangle$ and $|\phi\rangle$ as they amount to relative phases in the sum.
Given such abundance of possible paths, what would be the most natural choice of a path between $|\psi\rangle$ and $|\phi\rangle$? From a geometric perspective, the shortest (or geodesic) path along a great circle is a great choice.
If $\vec{r}, \vec{s} \in \mathbb R^3$ are the Bloch vectors of the two states, the line segment in $\mathbb R^3$ between them is described by $p \vec{r} + (1-p) \vec{s}$ where $p \in [0,1]$. Normalizing this Bloch vector to a unit vector gives you a geodesic path on the surface of the Bloch sphere.
To describe this in terms of density matrices, recall that the density matrix associated to Bloch vector $\vec{r} = (x,y,z)$ is given by $$ \rho(\vec{r}) = \frac{1}{2} (I + x X +y Y + z Z) = \frac{1}{2} \begin{pmatrix} 1+z & x-iy \\ x+iy & 1-z \end{pmatrix}. $$ You can check that $$ \operatorname{Tr}\rho(\vec{r})^2 = \frac{1}{2} (1 + \vec{r} \cdot \vec{r}) = \frac{1}{2} (1 + x^2 + y^2 + z^2), $$ meaning that $$ \|\vec{r}\| = \sqrt{x^2 + y^2 + z^2} = \sqrt{2\operatorname{Tr}\bigr(\rho(\vec{r})^2\bigl)-1}. $$ For any non-zero vector $\vec{r} \in \mathbb R^3$, the density matrix of a normalized pure state in direction $\vec{r}$ can be obtained as follows: $$ \rho\biggl(\frac{\vec{r}}{\|\vec{r}\|}\biggr) = \frac{1}{2} \Biggl( I + \frac{2\rho(\vec{r})-I} {\sqrt{2\operatorname{Tr}\bigr(\rho(\vec{r})^2\bigl)-1}} \Biggr). $$ Thus, given two arbitrary density matrices $\rho$ and $\sigma$ of pure qubit states, the geodesic between them is given by $$ \frac{1}{2} \Biggl( I + \frac{2(p\rho+(1-p)\sigma)-I} {\sqrt{2\operatorname{Tr}\bigl((p\rho+(1-p)\sigma)^2\bigr)-1}} \Biggr). $$ where $p \in [0,1]$. This gives you an unambiguous (while somewhat lengthy) way to interpolate between two pure qubit states.
You can also express this in terms of pure states. The key here is to note that the expression $$ a |\psi\rangle + b \frac{\langle\phi|\psi\rangle}{|\langle\phi|\psi\rangle|} |\phi\rangle $$ describes the same state regardless of the global phases of $|\psi\rangle$ and $|\phi\rangle$. By appropriately choosing $a$ and $b$ and normalizing the resulting state, you can also unambiguously interpolate between them.
Regarding your original question, the problem you are interested in is known as coherent state addition. You can look up Theorem 14 in my paper to see how it can be achieved, given enough copies of both states.
The main idea is to make one of the states evolve according to the Hamiltonian given by the commutator $i[|\phi\rangle\langle\phi|,|\psi\rangle\langle\psi|]$ of the two states. This can be achieved by repeatedly combining the two states using the swap test (this is reminiscent to what eqb proposed in their answer). The intuition for why this works is similar to Grover's algorithm -- you are trying to implement a rotation in the 2-dimensional subspace spanned by the two states and rotate one of them somewhere in between the two.
A counterexample that shows that this is not possible in the general case (here I am neglecting post-selection possibility discussed in the comments of the question and in the accepted answer):
$$ C_1 = X \qquad C_2 = -X$$
Or one can take $C_2 = R_y(- \pi)$ and all mentioned below equations will reamin true. So:
$$C_1 |0\rangle = |1\rangle = |\psi \rangle \qquad C_2 |0\rangle = -|1\rangle = |\phi \rangle$$
But we want to construct a circuit $C_3$:
$$C_3|0\rangle = \frac{|\psi \rangle + |\phi \rangle}{\sqrt{2}} = \frac{|1 \rangle - |1 \rangle}{\sqrt{2}} = 0$$
So $C_3$ is not unitary ($C_3$ doesn't preserve the length of the vector) and this shows that for general $C_1$ and $C_2$ this kind of circuit cannot be created.
Another counterexample without global phase ambiguity between $C_1$ and $C_2$:
$$C_1 = R_y(\frac{\pi}{3}) \qquad C_2 = R_y(-\frac{\pi}{3})$$
For this case;
$$C_1 |0\rangle = R_y(\frac{\pi}{3}) |0\rangle= \cos{\frac{\pi}{6}}|0\rangle + \sin{\frac{\pi}{6}}|1\rangle |\psi \rangle \\ C_2 |0\rangle = R_y(-\frac{\pi}{3}) |0\rangle = \cos{\frac{\pi}{6}}|0\rangle - \sin{\frac{\pi}{6}}|1\rangle = |\phi \rangle$$
We want to construct $C_3$:
$$C_3 |0\rangle = \frac{|\phi \rangle + |\psi\rangle}{\sqrt{2}} = \frac{\sqrt{3}}{\sqrt{2}} |0\rangle$$
And because $\frac{\sqrt{3}}{\sqrt{2}} \ne 1$, $C_3$ is not a unitary and cannot be constructed.
In some cases $C_3$ unitary exist: e.g. for $C_1 = X$ and $C_2 = Z$ we can have $C_3 = H = \frac{X + Z}{\sqrt{2}}$.
