Удалось найти фундаментальную работу по численным методам в гравитации. Собственно, я сам собирался сделать тоже самое

, но сначала решил погуглить по предполагаемым ключевым словам "разностный тензор кривизны", и вот, это оказалось уже сделано в 2015 году:
http://semr.math.nsc.ru/v12/p973-990.pdfЕ. Ю. Деревцов
Разностная аппроксимация ковариантной производной и других операторов и геометрических объектов, заданных в римановой области
Сибирские электронные математические известия
Siberian Electronic Mathematical Reports
http://semr.math.nsc.ruТом 12, стр. 973–990 (2015) УДК 517.9
DOI 10.17377/semi.2015.12.084
Поступила 3 августа 2015 г., опубликована 15 декабря 2015 г
Смысл работы в следующем.
В Римановой геометрии есть ряд тождеств:



А так же есть тождества Биянки и тождества Гильберта для тензора кривизны Римана.
При переходе к разностной схеме все тождества Римановой геометрии должны остаться тождествами. То есть для
разностной связности и для
разностного тензора кривизны разностные аналоги непрерывных тождеств Римановой геометрии должны удовлетворяться тождественно. Если вдруг какая-то взятая "с потолка" разностная схема для уравнений ОТО не уважает тождества Римановой геометрии, то она моделирует всё что угодно, но только не ОТО.
Оказывается, в разностном случае для ковариантного разностного дифференцирования ко- и контра- вариантных тензоров (формулы (1) и (2)) нужно использовать немножко разные связности (совпадающие в непрерывном случае).