I’m not surprised. There is challenge as we need to integrate densities and not log-densities. I did get some warnings about the accuracy already with 1D. One possibility is to evaluate first the log density with the two parameters set to their mean values, and then subtract that log-density from the log-densities inside integrand functions. Then the max density would be 1, and there is a smaller change of over/underflows. After getting the integration result and taking lof of that, the offset would be summed back.
You are integrating from -Inf to Inf, so I guess 1.79769e+308 is approximating Inf. You could try with a finite range, if you have a good guess for that (you can get the guess from normal distributions, maybe something like ± 9*sd)