The Galois group of $F/\mathbb Q$ is cyclic of order $6$ as you described above, and if $\sigma$ is a generator for this group, then since the unique subgroup of order $2$ is $\{1, \sigma^2, \sigma^4 \}$, we can assume without loss of generality that $\sigma^2 = u$. If you compute $\sigma$, you can compute $\sigma^2$ and explicit $E$ as the fixed field of $\sigma^2$.
Now the extension $F/\mathbb Q$ is simple, so $\sigma$ is characterized entirely by where it maps $f$ (let me write $f = \zeta_7$ for clarity). It suffices to find a primitive root of $(\mathbb Z/7\mathbb Z)^{\times}$ and we got it. After some testing you realize that $3^6 = 1$ but $3^2, 3^3 \neq 1 \pmod 7$, so $3$ is a primitive root. This means that $\sigma(\zeta_7) = \zeta_7^3$ is a possibility for $\sigma$ as a generator of the Galois group.
Now $\sigma^2(\zeta_7) = \zeta_7^9 = \zeta_7^2$ (as promised by $\sigma^2 = u$ and $u(f) = f^2$), so we're trying to find the fixed field corresponding to the group $\{ 1,\sigma^2,\sigma^4\}$.
Because of the simple relation $\zeta_7^7 = 1$, it is not hard to solve the equation $\sigma^2(x)-x = 0$ for $x \in F$ ; write
$$
x = a_0 + a_1 \zeta_7 + \cdots + a_5 \zeta_7^5
$$
and apply $\sigma^2$ to it ; the equation $\sigma^2(x) = x$ will give you linear conditions on the coefficients. I'm afraid finding generators involves linear algebra. If you have the patience of writing down the equation by hand, you'll see the computations do not require any smart-ideas and you get
$$
E = \mathbb Q(\zeta_7 + \zeta_7^2 + \zeta_7^4).
$$
The minimal polynomial can be computed, it is
$$
(X - (\zeta_7 + \zeta_7^2 + \zeta_7^4))(X - (\zeta_7^3 + \zeta_7^5 + \zeta_7^6)) = X^2 + X + 2.
$$
The roots of this polynomial are $\frac{-1 \pm \sqrt{-7}}2$, so you can write
$$
E = \mathbb Q(\sqrt{-7}).
$$
Note : The magical reason why I get only one canonical element (i.e. $\zeta_7 + \zeta_7^2 + \zeta_7^4$) as a generator for $E/\mathbb Q$ is because $[E:\mathbb Q] = 2$, so I know the extension can only be generated by one element and that a basis for $E/\mathbb Q$ will be given by $\{1,x\}$ for some $x \in E$. This means when I tried to compute the fixed field of $\sigma^2$, I know the term for $1$ would be irrelevant and the rest would only give me one term.
Interesting remark : You know that $\mathbb Z / 6 \mathbb Z$ has only one subgroup of index $2$ and that the extension $F/\mathbb Q$ is Galois, so by the Galois correspondence $F$ has only one subfield of order $2$ over $\mathbb Q$. The element $\zeta_7 + \sigma^2(\zeta_7) + \sigma^4(\zeta_7)$ is invariant by $\sigma^2$ and is not in $\mathbb Q$, so you could've skipped the linear algebra computations and know already that $E = \mathbb Q(\zeta_7 + \zeta_7^2 + \zeta_7^4)$. We don't escape the computation of the minimal polynomial to simplify this though. (These ideas of averaging are often used in invariant theory ; I averaged the orbit of $\zeta_7$ over the subgroup $\{1,\sigma^2,\sigma^4\}$ to find the invariant substructure I wanted.)
Hope that helps,
The thing is, when you prove the primitive element theorem, namely that $L:=F(\alpha,\beta)=F(\gamma)$ for some $\gamma\in L$, you need separability for only one of the roots, say $\beta$ (this is the only case that matters, for you can reduce inductively to it).
Take $\alpha_i, i=1,\dots, r$ and $\beta_j, j=1,\dots, s$ to be the distinct roots of minimal polynomials of $\alpha,\beta$ respectively in a splitting field. Since you can assume $F$ to be infinite (otherwise the proof is very easy), you can find $c\in F$ such that $\theta:=\alpha+c\beta$ differs from any $\alpha_i+c\beta_j$, the elements of form $\frac{\alpha-\alpha_i}{\beta-\beta_j}$ being finite.
If $\mu\in F[x]$ is the minimal polynomial of $\alpha$, you have that $\overline{\mu}(x):=\mu(\theta-cx)\in F(\theta)[x]$ verifies $\overline{\mu}(\beta)=0$ and $\overline{\mu}(\beta_j)\neq 0$ for all $\beta_j\neq\beta$. By the separability of $\beta$ on $F$ (and hence of $F(\theta)$), if $\nu(x)\in F[x]$ is the minimal polynomial of $\beta$, you obtain that $\text{gcd}(\nu(x),\overline{\mu}(x))= x-\beta$ in a splitting field $\overline{F}\supset F(\theta)$.
But $\text{gcd}$'s do not depend on the extension, and therefore you need $\beta\in F(\theta)$, which easily holds $\alpha\in F(\theta)$. Hence, $F(\theta)=F(\alpha,\beta)$.
Best Answer
Yes, the argument (essentially) works. Let me phrase it in a more careful way, without any assumption on the degree. Let $F_i:=E(\alpha_i)$ and view the $E_i$ as subfields of their Galois closure. There is an equivalence relation on the roots, namely, $\alpha_i \sim \alpha_j$ if $F_i = F_j$. The Galois group $G \subseteq S_n$ acts on the roots, and preserves this equivalence. In particular, either
No two roots generate the same field
All the roots generate the same field (and thus $F$ is Galois of degree $n$).
$G$ is an imprimitive subgroup of $S_n$ (i.e. preserves a non-trivial decomposition of $1,\ldots n$). This is clear, because the decomposition is given explicitly by the equivalence relation. In particular, this implies that $G$ is a subgroup of $S_A$ wreath $S_B$ for some $A,B \ne 1$ with $n = AB$.
If $n$ is prime, then, since $G$ acts transitively, 3 cannot occur, so "not 1" implies 2.