Cosmological hydrodynamic simulations of large scale structure in the universe have shown that accretion shocks and merger shocks form due to flow motions associated with the gravitational collapse of nonlinear structures. Estimated speed and curvature radius of these shocks could be as large as a few 1000 km/s and several Mpc, respectively. According to the diffusive shock acceleration theory, populations of cosmic-ray particles can be injected and accelerated to very high energy by astrophysical shocks in tenuous plasmas. In order to explore the cosmic ray acceleration at the cosmic shocks, we have performed nonlinear numerical simulations of cosmic ray (CR) modified shocks with the newly developed CRASH (Cosmic Ray Amr SHock) numerical code. We adopted the Bohm diffusion model for CRs, based on the hypothesis that strong Alfven waves are self-generated by streaming CRs. The shock formation simulation includes a plasma-physics-based 'injection' model that transfers a small proportion of the thermal proton flux through the shock into low energy CRs for acceleration there. We found that, for strong accretion shocks, CRs can absorb most of shock kinetic energy and the accretion shock speed is reduced up to 20%, compared to pure gas dynamic shocks. For merger shocks with small Mach numbers, however, the energy transfer to CRs is only about 10-20% with an associated CR particle fraction of 10-3. Nonlinear feedback due to the CR pressure is insignificant in the latter shocks. Although detailed results depend on models for the particle diffusion and injection, these calculations show that cosmic shocks in large scale structure could provide acceleration sites of extragalactic cosmic rays of the highest energy.