$\def\cX{\mathcal X}\def\cY{\mathcal Y}$The answer to the secondary question is an easy YES: for example, given an $L$-structure $\cX$ with domain $[n]$, let
$F(\cX)$ be the collection of pairs $(\cY,i)$ such that $\cY$ is an $L$-structure on $[n]$ isomorphic to $\cX$, $i\le n\log n$ is an integer, and the $i$th bit of the description of the lexicographically first isomorphism between $\cX$ and $\cY$ is $1$. Note that it is enough to have oracle access to *one* of $F(\cX)$ and $F(\cY)$; we do not need both.

If we want to make it depend less on details of the representation of the output, let e.g. $F(\cX)$ consist of pairs $(\cY,f)$ such that $\cY$ is an $L$-structure on $[n]$, 
$f$ is a partial self-map of $[n]$, and there exists an isomorphism $F\colon\cX\simeq\cY$ such that $F\supseteq f$. (It suffices to consider $f$ being the identity on an initial segment of the domain.) Using this, we can construct an isomorphism between $\cX$ and $\cY$ by fixing it one element at a time, querying, for a given $a\notin\operatorname{dom}(f)$, the existence of an isomorphism extending $f\cup\{(a,b)\}$ for each $b$ until we are successful.

We can express this by the truth of an existential first-order formula: if $f=\{(a_i,b_i):i<k\}$, let $\delta(\vec b,\vec y)$ be the (conjunction of the) atomic diagram of $(\cY,\vec b)$; then $f$ extends to an isomorphism $\cX\simeq\cY$ iff $\mathcal X\models\exists\vec y\,\delta(\vec a,\vec y)$. Thus, you can also take for $F(\mathcal X)$ a truth predicate for existential first-order formulas with parameters from $\cX$.

The same kind of argument also works for the main question: given structure $\cX$ and $\cY$ with domain $[n]$, and a partial self-map $f$ on $[n]$, let $\delta_\cX(\vec x)$ be the atomic diagram of $\mathcal X$ in variables $\{x_i:i\in[n]\}$, and $\delta_\cY(\vec y)$ be the atomic diagram of $\cY$ in variables $\{y_i:i\in[n]\}$. Then $f$ extends to an isomorphism $\cX\simeq\cY$ iff the existential second-order sentence (with no extra-logical symbols)
$$\exists F\,\exists\vec R\,\exists\vec x\,\Bigl(\delta_\cX(\vec x)\land\delta_\cY(F(\vec x))\land\bigwedge_{f(i)=j}F(x_i)=x_j\Bigr)$$
is true in some/every structure with $n$ elements, where $\vec R$ quantifies over all relations and functions of $L$ as used in $\delta_\cX$ and $\delta_\cY$, and $F(\vec x)$ denotes the sequence of values $F(x_0),F(x_1),\dots$. Again, using this oracle, we can construct an isomorphism $\cX\simeq\cY$ one element at a time; we only use the oracle for one of $\mathrm{Th}_2(\cX)$ or $\mathrm{Th}_2(\cY)$ rather than their join.

The only issue is whether this can be done in *linear* time. Since this may depend on all kinds of details, let me assume that we are using the standard multi-tape Turing machine model with a read-only input tape, write-only output tape, several work tapes, and an oracle query tape. The content of the oracle query tape is not erased or otherwise modified after making an oracle query.

We can write down the $\exists F\,\exists\vec R\,\exists\vec x\,\bigl(\delta_\cX(\vec x)\land\delta_\cY(F(\vec x))\land\dots$ part of the formula on the query tape in time linear in the size of $\mathcal X$ and $\mathcal Y$ (as the diagrams $\delta_\cX$ and $\delta_\cY$ are basically just copies of the tables of the relations and functions). Then we write down $\dots F(x_0)=x_0\bigr)$ and query the oracle; if the answer is negative, we modify it to $\dots F(x_0)=x_1\bigr)$, and so on, until we get a positive answer. Each such change is constant time if the variable indices are written in unary; if they are written in binary, it’s amortized constant time. Anyway, it takes time $O(n)$ to find the image of $0$; we copy it to the output tape, and go on to add $\dots\land F(x_1)=x_0$ on the query tape, and continue. In this way, we construct the isomorphism in time $O(n^2)$ on top of the initialization part which was linear in the size of the input.

Thus, if $L$ contains at least one at least binary relation or function symbol, the algorithm as a whole works in linear time, as the size of the input is $\Omega(n^2)$.

What if all symbols in $L$ are unary? Well, first of all, if we have only unary relations and constants, then the size of the input is $O(n)$, whereas the size of the output is $\Omega(n\log n)$. Thus, the task is impossible to do in linear time with whatever oracle, as a linear-time machine does not even have the time to write down the result. (So the claim in the first paragraph of the question is wrong, by the way.) The best we can hope for in this case is to make it work in time $O(n\log n)$; that is, to have an algorithm that works in time linear in the combined size of the input and output.

I’m not sure whether that is possible to arrange. We can reduce the number of oracle queries to $O(n\log n)$ if we use binary search to find the image of $i\in[n]$: i.e., instead of trying whether the image can be $0$, $1$, $2$, ..., we first determine the first bit of the image, then the second bit, and so on. To make sure that the individual queries do not take too much time to construct, we can introduce predicates for elements with a given bit set: the formula will now start with
$$\exists F\,\exists\vec R\,\exists\vec x\,\exists P_0,\dots,P_{\log n}\,\Bigl(\delta_\cX(\vec x)\land\delta_\cY(F(\vec x))\land\bigwedge_{\substack{i\in[n]\\j\le\log n}}(P_j(x_i))^{b(i,j)}\land\dots,$$
where $b(i,j)\in\{0,1\}$ is the $j$th bit of $i$, and $\phi^1=\phi$, $\phi^0=\neg\phi$. But just writing down the indices of all the new predicates in binary in these $n\log n$ conjuncts takes $O(n\log n\log\log n)$ bits. One can then arrange the whole computation to take time $O(n\log n\log\log n)$ (plus linear in the size of the input in case of nonunary symbols), but I’m not sure how to reduce it all the way down to $O(n\log n)$.