The crystal plasticity finite element model (CPFEM) is a significant tool in the integrated computational materials engineering (ICME) toolboxes that bridges between microstructures and materials properties relationship. However, to establish the predictive capability, one needs to calibrate the underlying constitutive model, verify the numerical solution, and validate the model prediction against experimental data. Bayesian optimization (BO) has stood out as a gradient-free efficient global optimization algorithm that is capable of calibrating constitutive models for CPFEM. To cope with the stochastic nature of representative volume elements in microstructure, the Monte Carlo estimator is employed to approximate the loss function over a microstructure ensemble. In this work, we propose a high-throughput asynchronous parallel workflow to calibrate constitutive models through Bayesian optimization, where multiple CPFEM simulations are performed concurrently. The proposed method is demonstrated for stainless steel 304L, tantalum, and Cantor high-entropy alloy.