General Relativity and Computing with Tensors

How complicated or expensive the problem/computations become at the end highly depends on what you do and what one really is interested in.

In the scope of your question lets look at a simple problem (which involves numerical solutions) and how one would really tackle it: A spherical symmetric, static, cold compact star composed of an ideal fluid. This is the setting of the Tolman-Oppenheimer-Volkoff (TOV) equation. After some algebra, which works in the symmetries one is left with the task of solving 3 differential equations: two for the only two open metric potentials and 1 for the pressure gradient. The Einstein equations come down to just 3 differential equations. I can bring them to this state just by working in my symmetries. One does not need any specific values of the tensor fields. To bring down/simplify expressions one does not need the functional forms: eg. I can lower or raise indices with the metric where and when ever I want for that I do not need to know how my metric explicitly looks like in the end: $x_mu=g_{munu}x^nu$ holds for all $x_muin mathbb{R}^4$ and I can write this down without having any clue about the metric or how the metric potentials look like. All preliminary tensor algebra looks works like that.

So to completely describe the star and its spacetime one only needs to solve those 3 differential equations numerically. What one is left with after are three functions $P(r)$, $nu(r)$ and $lambda(r)$. $P(r)$ describes the pressure gradient of the star and with an equation of state this specifies the configuration completely and $nu(r)$ and $lambda(r)$ are the metric potentials from the line element: $$ds^2=-e^{nu(r)}dt^2+e^{lambda(r)}dr^2+r^2(dtheta^2+sin^2theta dphi^2).$$

Now that we know all missing parts we can compute anything related to the problem: geodesics, total internal energy, surface gravity everything one can think of. The total internal energy would be something like $$M_{TW}=int (T_{i}^{~i}-T_{t}^{~t})sqrt{-g}dV =4pi int_0^R dr~ r^2e^{(nu(r)+lambda(r))/2}(rho(r)+3P(r)).$$ Computing that numerically is straight forward once we have our metric potentials and pressure gradient as numerical solutions. So that we do not know analytical solutions for our metric potentials and pressure gradients is not a problem in the end because we really do not need to have them analytically.

If one is really interested in high end numerical General relativistic calculations in 4D I would recommend literature on Numerical Relativity. Or you could look at software packages/codes like LORENE or Whisky. Those codes implement tensor fields and operations on and with them. Those codes are able to solve the Einstein equations (+Maxwell,...) on spacetime grids numerically exact. Depending on the problem at hand however one will always try to work in the symmetries to reduce the dimensionality as much as possible. In the scope of those codes the Einstein equations in general do not decouple and very involved numerical methods are necessary to solve those coupled differential equations.

Thinking about it in abstract geometrical terms might help. So, you have a (pseudo-)Riemannian manifold but you don't have to consider coordinates.

In calculations, you're considering a point on the manifold, let's say $mathbf{x}$. All quantities are either evaluated at that point or they are integrals over specified regions. One cannot trivially compare values of fields at different points and that's why we use Christoffel symbols $mathbf{Gamma}$, i.e. connections. They allow us to take into account the metric structure of the manifold and go between different coordinate patches.

Raising, lowering and contracting indices is accomplished by multiplication with $g^{mu nu} (x)$. So what happens there?

Consider a vector-valued function $A^mu (x)$. This function takes a point on the manifold, parametrized by $x$, and it maps it to a contravariant vector at that point, i.e. to an element of the cotangent vector space at that point.

Similarly, $A_{mu} (x)$ maps the point to a covariant vector, i.e. an element of the tangent vector space at that point.

The metric tensor, on the other hand, is perhaps best viewed as a function that takes two arguments and spits out a scalar. Specifically, it takes a tangent vector and a cotangent vector and maps them to a scalar. The thing is, you have a different (co)tagent space at each point on the manifold (they constitute the (co)tangent bundle, by the way) and you obviously can't work with vectors from different vector spaces. That's why $mathbf{g}$ is a function of position, you have to specify a point and both vectors have to be defined in spaces at that point on the manifold.

With that being said, raising and lowering simply amounts to looking at the metric with one slot empty, i.e. $$g^{mu nu}A_nu leftrightarrow mathbf{g} ( mathbf{A}, , , )$$ $$g_{mu nu}A^nu leftrightarrow mathbf{g} ( , , , mathbf{A} )$$

Contraction should be obvious at this point, although it's just slightly more complicated to talk about it.

The point (no pun intended) is, those operations are all taken at the same point.