Does $R^2$ increase with the number of variables no matter what model we are adjusting? I understand this is true for OLS, but i am not sure if this is true for every other model, like ridge or a NN.

The question is a bit vaguely formulated (how is one increasing the number of variables?, what sort of loss function is your model trying to minimize?). In the following result, I try to state a set of precise conditions under which the answer to your question is "yes". I then discuss why my result is true under the stated conditions as well as when my stated conditions fail to give you a sense for what sort of "pathologies" as alluded to in @whuber's comment one would need for $R^2$ to not increase.
Claim: Let $X_1$ be a vector of predictors and let $X_2$ be a vector of predictors containing $X_1$ as a sub-vector. Let $Y$ be the dependent variable whose value is being predicted. Let $\theta_1\in \mathcal H_1$ be the parameters of a model using $X_1$ as predictors such that when the model parameter is $\theta_1$, our predictions of $Y$ given $X_1$ are $m_1(X_1;\theta_1)$ and let $\theta_2 \in \mathcal H = \mathcal H_1 \cup \mathcal H_2$ be the parameters of a model using $X_2$ as predictors such that when the model parameter is $\theta_2$, our prediction of $Y$ given $X_2$ is $m_2(X_2,\theta_2)$. Finally, assume that our models are nested in the sense that if $\theta_2\in\mathcal H_1$, then $m_1(X_1,\theta_2) = m_2(X_2,\theta_2)$ whenever $X_1$ has the same values as $X_2$ in all shared commponents. Finally, suppose that our "best" predictors are chosen via least squares, i.e. for $i=1,2$,
$$\theta_i^* = \mathrm{argmin}_{\theta_i}\, \mathbb E[(Y-m_i(X_i,\theta_i))^2]$$
Define the population $R^2$ of model $i$ by the quantity
$$R^2_i \equiv 1 - \frac{\mathbb E[(Y-m_i(X_i,\theta_i^*))^2]}{\mathrm{Var}(Y)}$$
Then $R_1^2 < R_2^2$.
The claim required a lot of setup to state, but the proof is simple. Note that when one is fitting a least-squares model, the objective being minimized is exactly the numerator of the object being subtracted in $R^2$, and the denominator depends only on properties of $Y$. Moreover, when fitting using predictors $X_2$, our assumption that the models are nested implies that a valid possible value of $\theta_2$ is one that effectively ignores the extra information in $X_2$ and predicts according to $\theta_1^*$. Thus, the fact that $\theta_2^*$ minimizes the objective implies that one can only do better with a richer set of predictors. Formally, by definition of minimization, we have
$$\mathbb E[(Y-m_2(X_2,\theta_2))^2] \leq \mathbb E[(Y-m_2(X_2,\theta_1^*))^2] = \mathbb E[(Y-m_1(X_1,\theta_1^*))^2]$$
and the rest follows from some simple algebra.
So now, what are the possible violations to my long list of assumptions. First, it is possible that $X_2$ has a larger number of variables than $X_1$, but there are some variables in $X_1$ that are highly predictive but not in $X_2$. In this case, model 2 might have a lower population $R^2$. Second, it is possible that the objective function used to choose the "best" $\theta_i$ is not the least-squares objective. This is essentially @Dave's example, and there again, there is no reason in general to assume that a $\theta_i$ chosen to fit some other objective happens to also fit the least-squares objective well too.
