So we try to find an approximation to $p(θ|x)$ by variational or
sampling methods.
This may need to be clarified a little. Variational methods approximate $p(
\theta |x)$ by the closest possible approximation within some class of distributions. Sampling methods do not really set out to approximate $p(\theta|x)$ (although one could do that). Instead, usually one is interested in expectations or quantiles of $p(\theta | x)$, and they set out to estimate that. The key here is the word estimate, instead of approximate. You can find a more elaborate difference between the two methods in my answer here.
1) In order to do sampling we need to calculate $p(x|\theta)p(\theta)$
for a specific value of $x$ and $\theta$, which gives us an
unnormalized scalar value, and we know the exact shapes of
$p(x|\theta)$ and $p(\theta)$ distributions separately. Right?
Yes, usually in most cases, we have known likelihood $p(x|\theta)$ and a chosen known prior $p(\theta)$. To begin with, we assume that we can evaluate $p(x|\theta)p(\theta)$ in closed form for a given value of $x$ and $\theta$. This can sometimes be computationally expensive, if evaluating the likelihood is expensive.
2) Is it possible that we also do not know or cannot compute the
likelihood and prior? Can Approximate Bayesian Computation or other
likelihood-free inference methods be used for this?
Absolutely, it is possible that the likelihood specially cannot be evaluated (I haven't come across examples where the prior cannot be evaluated). In such cases, as you say, one can use Approximate Bayesian Computation, when there is no known distributional form of the likelihood, and data comes from a complex machine (for example climate models).
Alternatively, sometimes likelihoods themselves are intractable. This can happen for example in generalized linear mixed models, in problems with missing data etc. In such cases one can use pseudo-marginal MCMC methods, where inside the MCMC step, an unbiased estimate of the likelihood replaces the true likelihood evaluation.
3) In rejection sampling the proposal density has to cover the
posterior, right? How can we know if it does or not? We don't know the
exact shape since we don't know the normalization term.
Yes, in rejection sampling, the support of the proposed density must contain the support of the target density. This usually easy to determine, and is governed by the prior on $\theta$. For example if I assume the prior $\theta \sim Beta(\alpha, \beta)$, then I know that the support of $\theta$ is $(0,1)$. If my prior is $\theta \sim N(0,1)$, then I know that the support of $\theta$ is the full real line.
4) In MCMC this is not a problem since the normalization terms cancel
when dividing two posterior values (which are scalar). Therefore just
evaluating $p(x|\theta)p(\theta)$ is enough though we don't know the
exact shape of the whole distribution. Is this correct?
Yes.
5) If we sample the posterior, then we have $n$ samples $\theta_1, ..., \theta_n$. But we still don't have the distribution
$p(\theta|x)$, so how can we find this?
Like I wrote in the beginning, it depends on what quantities you are interested in. If you are interested in the posterior mean of $\theta$, then you are interested in
$$E_{\theta|x}[\theta] = \int \theta p(\theta|x)\,. $$
This is estimated by
$$\dfrac{1}{n} \sum_{t=1}^{n} \theta_t\,. $$
This is why MCMC is called as such (Markov chain Monte Carlo). The sampling is done using a Markov chain, and the estimation is done using Monte Carlo.
Similarly, all expectation, and quantiles can be estimated using such Monte Carlo estimators. There are also methods that estimate the normalizing constant using the samples, to give you an estimate of $p(\theta |x)$, like the paper here. However, in most cases, one is only interested in attributes of $p(\theta|x)$, which can directly be estimated using the MCMC samples.
6) For new data points we need to compute the posterior predictive
$p(x_{new}|D)=\int p(x_{new}|\theta,D)p(\theta|D)d\theta$ where
$D=\lbrace x_1,x_2,...,x_m \rbrace$ is the training set. How can we do
this using the posterior samples? Do we again need sampling to
calculate this integration?
Notice that the posterior predictive is nothing but an expectation with respect to the posterior distribution. In addition, note that given $\theta$, the new data point is independent of the old data. That is, $p(x_{new}| \theta, D) = p(x_{new} | \theta)$. So to estimate the posterior predictive
$$\dfrac{1}{n} \sum_{t=1}^{n} p(x_{new}|\theta_t,D) = \dfrac{1}{n} \sum_{t=1}^{n} p(x_{new}|\theta_t)\,, $$
as long as you can evaluate the likelihood.
7) What does computing and evaluating a distribution mean? What is the
difference?
Note first that there is a difference between distribution and the density of the distribution (read the discussion here. Evaluating a distribution, usually refers to evaluating the density at a particular point. For example, to use MCMC, you need to evaluate $p(\theta)$ at the proposed value of $\theta = \theta_{prop}$.
Computing a distribution would refer to finding the distribution, as in knowing the full density or cumulative density function of the distribution. So for example, while implementing MCMC, in every step we evaluate the posterior distribution (up to a normalizing constant), but at the end of the MCMC, if we can estimate the normalizing constant, we would have computed the distribution.
