hevgyrt / ocean_wave_tracing Goto Github PK

import numpy as np
import maplotlib.pyplot as plt
from ocean_wave_tracing.ocean_wave_tracing import Wave_tracing
 
# Defining some properties of the medium
nx = 100; ny = 100 # number of grid points in x- and y-direction
x = np.linspace(0,2000,nx) # size x-domain [m]
y = np.linspace(0,3500,ny) # size y-domain [m]
T = 250 # simulation time [s]
U=np.zeros((nx,ny))
U[nx//2:,:]=1
 
# Define a wave tracing object
wt = Wave_tracing(U=U,V=np.zeros((ny,nx)),
                       nx=nx, ny=ny, nt=150,T=T,
                       dx=x[1]-x[0],dy=y[1]-y[0],
                       nb_wave_rays=20,
                       domain_X0=x[0], domain_XN=x[-1],
                       domain_Y0=y[0], domain_YN=y[-1],
                       )
 
# Set initial conditions
wt.set_initial_condition(wave_period=10,
                              theta0=np.pi/8)
# Solve
wt.solve()

# Plot
fig, ax = plt.subplots();
pc=ax.pcolormesh(wt.x,wt.y,wt.U.isel(time=0));
fig.colorbar(pc)

for ray_id in range(wt.nb_wave_rays):
    ax.plot(wt.ray_x[ray_id,:],wt.ray_y[ray_id,:],'-k')

plt.show()

        vo = np.full(shape, np.nan, dtype=block.dtype)
        vo[:, inbounds] = block.values[:, ty.ravel(), tx.ravel()]

        vo = xr.DataArray(vo,
                          [
                              ("time", var.time.data),
                              ("latitude", target.y),
                              ("longitude", target.x),
                          ],
                          attrs=var.attrs,
                          name=var.name)

        # Positions in source grid
        # vo.attrs['x'] = target_x
        # vo.attrs['y'] = target_y

        vo.latitude.attrs['units'] = 'degrees_north'
        vo.latitude.attrs['standard_name'] = 'latitude'
        vo.latitude.attrs['long_name'] = 'latitude'

        vo.longitude.attrs['units'] = 'degrees_east'
        vo.longitude.attrs['standard_name'] = 'longitude'
        vo.longitude.attrs['long_name'] = 'longitude'
        vo.attrs['grid_mapping'] = target.proj_name
        vo.attrs['source'] = self.url
        vo.attrs['source_name'] = self.name

    ds = xr.Dataset()

    for var in sources.scalar_variables:
        logger.info(f'Searching for variable {var}')
        (d, v) = sources.find_dataset_for_var(var)

        if v is not None:
            logger.info(f'Extracting {var} from {d}')

            # Acquire variables on target grid
            vo = d.regrid(v, target, t0, t1)
            ds[vo.name] = vo
        else:
            logger.error(f'No dataset found for variable {var}.')

    for vvar in sources.vector_variables:
        varx = vvar[0]
        vary = vvar[1]
        logger.info(f'Searching for variable {varx},{vary}')
        (d, vx, vy) = sources.find_dataset_for_var_pair(varx, vary)

        if vx is not None:
            logger.info(f'Extracting {varx} and {vary} from {d}')

            # Acquire variables on target grid
            vox = d.regrid(vx, target, t0, t1)
            voy = d.regrid(vy, target, t0, t1)

            vox, voy = d.rotate_vectors(vox, voy, target)

            ds[vox.name] = vox
            ds[voy.name] = voy
        else:
            logger.error(f'No dataset found for variable {varx},{vary}.')