Taking a derivative, then squaring a two dimensional numerical function

Question

Context

I am interested in first finding an interpolating function of the solution to the linearly damped wave equation. Here is the solution to the LDWE with smooth square inital function :

ub[epsilon_] := 
Interpolation[{{{-5.944 - epsilon}, 0, 0}, {{-5.944 + epsilon}, 1,0}, 
{{5.944 - epsilon},     1, 0}, {{5.944 + epsilon}, 0, 0}}, 
"ExtrapolationHandler" -> {0 &, "WarningMessage" -> False}]

square = ub[0.0001]; initialWavePulseSquare[x_] := square[x];

exactSolSquare[x_, t_, c_, τ_] := {pEvaluate = Sqrt[1 - (x - s)^2/(c t)^2];
 limitInt1 = x - c t; limitInt2 = x + c t; 
 1/2 (initialWavePulseSquare[x - c t] + initialWavePulseSquare[x + c t])*Exp[-t/(2 τ)]  
 + 1/(4 c τ) Exp[-t/(2 τ)]*
 Quiet[NIntegrate[
  initialWavePulseSquare[s]*(BesselI[0, pEvaluate*t/(2 τ)] + 
  (1/pEvaluate)* BesselI[1, pEvaluate*t/(2 τ)] ) , {s, limitInt1, 
   limitInt2}]]}

 solutionRewritten[x_, t_] := exactSolSquare[x, t, 1, 25]

Here, I used

square = ub[0.0001]; 
initialWavePulseSquare[x_] := square[x] 

to define a "smooth" rectangular initial wave pulse. Then, I defined the solution to the LDWE as exactSolSquare[x, t, c, τ] where c and tau are speed and damping constant, respectively. Then, I was able to create a numerical two dimensional solution to the LDWE defined as solutionRewritten[x_, t_]

Objective

I want to ultimately take the integral of (Derivative[0,1][solutionRewritten][x,t])^2 from t=0 to t=300. However, I have hard time even getting the InterpolatingFunction.

Question

How can I create an interpolating function of two dimensional numerical data and evaluate the integral of (Derivative[0,1][solutionRewritten][x,t])^2 from t=0 to t=300.

Attempt 1

I tried something like this and MM threw a handful of errors:

Interpolation[Table[Evaluate[{{x, t}, solutionRewritten[x, t]}], {x, -200, 200, .5},    
{t, 0.01, 300, 0.5}] ]

I previously asked similar question before. But this time, I know the solution of NDSolve explicitly. Here is the link: Integrating Squared of Interpolating Function with respect to one variable

Attempt 2

I tried using a smaller boundary to see if the code even runs. And as expected, it does run and gives me an InterpolatingFunction.

solnInterpolation = ListInterpolation[Table[First[solutionRewritten[x, t]], {x, -20, 20}, {t, 0, 20}]]

However, it outputs the following error: Power::infy: "Infinite expression 1/0 encountered". I tried to suppress the error messages and re-evaluated the code with larger boundary as follows:

Off[Power::infy]; solnInterpolation = ListInterpolation[Table[First[solutionRewritten[x, t]], {x, -200, 200}, {t, 0, 300}]]

However, MM evaluates this line of code runs indefinitely and doesn't output an InterpolatingFunction. Is there a way to speed this process with pretty good precision?

Thank You very much for your help!