64 bit floats

wecopttool/core.py

@@ -725,15 +725,21 @@ def solve(self,
             scale_x_opt = scale_dofs([scale_x_opt], nstate_opt)
 
         # composite scaling vector
-        scale = jnp.concatenate([jnp.array(scale_x_wec), jnp.array(scale_x_opt)])
+        scale = jnp.concatenate([
+            jnp.array(scale_x_wec, dtype=jnp.float64),
+            jnp.array(scale_x_opt, dtype=jnp.float64),
+        ])
 
         # decision variable initial guess
         key = jax.random.PRNGKey(0) # could add key as input to select same initial guesses?
         if x_wec_0 is None:
-            x_wec_0 = jax.random.normal(key, [self.nstate_wec], dtype=np.float64)
+            x_wec_0 = jax.random.normal(key, [self.nstate_wec])
         if x_opt_0 is None:
-            x_opt_0 = jax.random.normal(key, [nstate_opt], dtype=np.float64)
-        x0 = jnp.concatenate([jnp.array(x_wec_0), jnp.array(x_opt_0)])*scale
+            x_opt_0 = jax.random.normal(key, [nstate_opt])
+        x0 = jnp.concatenate([
+            jnp.array(x_wec_0, dtype=jnp.float64),
+            jnp.array(x_opt_0, dtype=jnp.float64),
+        ]) * scale
 
         # bounds
         if (bounds_wec is None) and (bounds_opt is None):
@@ -2332,7 +2338,8 @@ def wave_excitation(exc_coeff: DataArray, wave: DataArray) -> ndarray:
             f"\n Wave direction(s): {(np.rad2deg(dir_w))} (deg)" +
             f"\n BEM direction(s): {np.rad2deg(dir_e)} (deg).")
 
-    return jnp.sum(wave_elev_fd*exc_coeff[:, sub_ind, :], axis=1)
+    result = wave_elev_fd.astype(jnp.complex128) * exc_coeff[:, sub_ind, :].astype(jnp.complex128)
+    return jnp.sum(result, axis=1)
 
 
 def hydrodynamic_impedance(hydro_data: Dataset) -> Dataset:

wecopttool/pto.py

@@ -235,7 +235,8 @@ def _fkinematics(self,
         assert f_wec_td.shape == (wec.nt*nsubsteps, wec.ndof)
         f_wec_td = jnp.expand_dims(jnp.transpose(f_wec_td), axis=0)
         kinematics_mat = self.kinematics(wec, x_wec, x_opt, wave, nsubsteps)
-        return jnp.transpose(jnp.sum(kinematics_mat*f_wec_td, axis=1))
+        result = kinematics_mat.astype(jnp.float64) * f_wec_td.astype(jnp.float64)
+        return jnp.transpose(jnp.sum(result, axis=1))
 
     def position(self,
         wec: TWEC,
@@ -362,7 +363,9 @@ def force_on_wec(self,
         assert force_td.shape == (1, self.ndof, wec.nt*nsubsteps)
         kinematics_mat = self.kinematics(wec, x_wec, x_opt, wave, nsubsteps)
         kinematics_mat = jnp.transpose(kinematics_mat, (1,0,2))
-        return jnp.transpose(jnp.sum(kinematics_mat*force_td, axis=1))
+        # Use float64 for higher precision in the multiplication and sum
+        result = kinematics_mat.astype(jnp.float64) * force_td.astype(jnp.float64)
+        return jnp.transpose(jnp.sum(result, axis=1))
 
     def mechanical_power(self,
         wec: TWEC,
@@ -425,7 +428,8 @@ def mechanical_energy(self,
             length.
         """
         power_td = self.mechanical_power(wec, x_wec, x_opt, wave, nsubsteps)
-        return jnp.sum(power_td) * wec.dt/nsubsteps
+        # Use float64 for higher precision in the sum
+        return jnp.sum(power_td.astype(jnp.float64)) * wec.dt/nsubsteps
 
     def mechanical_average_power(self,
         wec: TWEC,
@@ -572,7 +576,8 @@ def energy(self,
             length.
         """
         power_td = self.power(wec, x_wec, x_opt, wave, nsubsteps)
-        return jnp.sum(power_td) * wec.dt/nsubsteps
+        # Use float64 for higher precision in the sum
+        return jnp.sum(power_td.astype(jnp.float64)) * wec.dt/nsubsteps
 
     def average_power(self,
         wec: TWEC,

wecopttool/utilities.py

@@ -252,9 +252,9 @@ def calculate_power_flows(wec,
         P_e.append((1/4)*(Fe_FD_t@np.conj(U_FD) + U_FD_t@np.conj(Fe_FD)))
 
     power_flows = {
-        'Optimal Excitation' : -2* np.sum(np.real(P_max)),#eq 6.68
-        'Radiated': -1*np.sum(np.real(P_r)),
-        'Actual Excitation': -1*np.sum(np.real(P_e)),
+        'Optimal Excitation' : -2* np.sum(np.real(P_max), dtype=np.float64),#eq 6.68
+        'Radiated': -1*np.sum(np.real(P_r), dtype=np.float64),
+        'Actual Excitation': -1*np.sum(np.real(P_e), dtype=np.float64),
         'Electrical (solver)': P_elec,
         'Mechanical (solver)': P_mech,
                   }