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.

There are less common estimators that we can use and might want to use under some circumstances. For instance, if we estimate the parameters of a linear model by minimizing something other than square loss (sum of squared residuals), we might wind up with the square loss increasing as we add parameters, since the optimization did not know to try to minimize square loss.
In the simulation below, I have an estimator that minimizes a quantile loss to estimate the coefficients. The result is the smaller model has a smaller sum of squared residuals and, consequently, a larger $R^2$.
library(quantreg)
set.seed(2022)
N <- 1000
x1 <- rnorm(N)
x2 <- rnorm(N)
yh <- x1
e <- rnorm(N, 0, 1)
y <- yh + e
L <- lm(y ~ x1 + x2)
X1 <- cbind(x1)
X2 <- cbind(x1, x2)

quantile_r2 <- function(X, y){
  
  L <- quantreg::rq(y ~ X, tau = 0.6)
  
  preds <- predict(L)
  
  SSTotal <- sum(
    (
      y - mean(y)
    )^2
  )
  
  SSResiduals <- sum(
    (
      y - preds
    )^2
  )
  
  R2 <- 1 - SSResiduals/SSTotal
  
  return(R2)
}
quantile_r2(X1, y) - quantile_r2(X2, y)

$$
R^2 = 1-\dfrac{
\underset{i=1}{\overset{n}{\sum}}\left(
y_i - \hat y_i
\right)^2
}{
\underset{i=1}{\overset{n}{\sum}}\left(
y_i - \bar y
\right)^2
}
$$
Overall, it depends on what exactly you're doing if you can claim that more parameters necessarily results in higher in-sample $R^2$.
