def loss_ic(model, x_ic, rho_ic, u_ic, p_ic):
    """
       Loss function for initial condition.
    """
    
    y_ic = model.net(x_ic)  # Initial condition
    
    rho_ic_nn = y_ic[:, 0]
    p_ic_nn = y_ic[:, 1]
    u_ic_nn = y_ic[:, 2]

    # Loss function for the initial condition
    loss_ics = ((u_ic_nn - u_ic) ** 2).mean() + \
           ((rho_ic_nn - rho_ic) ** 2).mean()  + \
           ((p_ic_nn - p_ic) ** 2).mean()

    return loss_ics