I want to define two different regions with different properties, however, the following error is coming and I don't know how can I sort it out

Question

import numpy as np
import ComPASS
from ComPASS.utils.units import *
from ComPASS.utils.grid import grid_center
from scipy.interpolate import griddata

#Set Output information
ComPASS.set_output_directory_and_logfile(__file__)

# Parameters
pres = 30. * MPa            # initial reservoir pressure
Tres = degC2K( 200. )        # initial reservoir temperature - convert Celsius degrees to Kelvin degrees

#Tinjection = degC2K( 70. )  # injection temperature - convert Celsius to Kelvin degrees

#Qm = 1200. * ton / hour      # production flowrate

#interwell_distance = 1.2 * km # distance between wells

#n_COX = 0.15
#n_CCT = 0.2      # reservoir porosity

#k_COX= 1E-18 
#k_CCT = 1.5E-12         # reservoir permeability in m^2

#K_reservoir= 1.5
K_reservoir = 2             # bulk thermal conductivity in W/m/K

g = 9.81                    #specific gravity in m/s^2

#UCG_heat_flux = 2.0  # W/m2
#rho_rock = 2000 # kg/m^3 rock specific mass
#cp_rock = 800  # J/K/kg specific heat capacity

#Meshes and grids
Lx, Ly, Lz = 4000.0, 2500.0, 1500.0
Ox, Oy, Oz = -1500.0, -1000.0, -500.0
nx, ny, nz = 50, 30, 20
#
grid = ComPASS.Grid(
    shape=(nx, ny, nz),
    extent=(Lx, Ly, Lz),
    origin=(Ox, Oy, Oz),
)

print('grid shape: ',grid.shape)

# Define the rock types
def define_rock_mass_boundaries():
    cell_centers = simulation.compute_global_cell_centers()

    # Determine the x-coordinate ranges for COX and CCT regions
    min_x_cox = 0
    max_x_cox = 2499
    min_x_cct = 2500
    max_x_cct = 4000


    # Identify cells within the COX and CCT regions based on x-coordinate
    cox_indices = np.where((cell_centers[:, 0] >= min_x_cox) & (cell_centers[:, 0] <= max_x_cox))
    cct_indices = np.where((cell_centers[:, 0] >= min_x_cct) & (cell_centers[:, 0] <= max_x_cct))
    
    print ("cell_centers shape:",cell_centers.shape)
    print ("cct indices:",cct_indices)
    print ("cox indices:",cox_indices)
    
    return [cox_indices,cct_indices]


def porosity(cox_indices,cct_indices):
    porosity_grid = np.ones(cell_centers.shape)  
    porosity_coefficient_cox = 0.15  # Adjust the porosity value for COX region as desired
    porosity_coefficient_cct = 0.3  # Adjust the porosity value for CCT region as desired
    porosity_grid[cox_indices] = porosity_coefficient_cox
    porosity_grid[cct_indices] = porosity_coefficient_cct

    return porosity_grid
    
def permeability(cox_indices,cct_indices):    
    permeability_grid = np.ones(cell_centers.shape)  
    perm_coefficient_cox = 1e-14 
    perm_coefficient_cct = 1e-18 
    permeability_grid[cox_indices] = perm_coefficient_cox  # Use the desired permeability value for COX region   
    permeability_grid[cct_indices] = perm_coefficient_cct  # Use the desired permeability value for CCT region
    return permeability_grid
   
def density(cox_indices,cct_indices):    
    density_grid = np.ones(cell_centers.shape)  # Assuming initial density of 1 (no variation)
    density_coefficient_cox = 2.65 
    density_grid[cox_indices] = density_coefficient_cox
    density_coefficient_cct = 2.4
    density_grid[cct_indices] = density_coefficient_cct
    return density_grid
    
def thconductivity(cox_indices,cct_indices):
    thconductivity_grid=np.ones(cell_centers.shape)   
    thconductivity_cox=2.1
    thconductivity_cct=4
    thconductivity_grid[cox_indices] = thconductivity_cox
    thconductivity_grid[cct_indices] = thconductivity_cct

    return thconductivity_grid

#Load the physics of the system
#Definition
simulation = ComPASS.load_physics("immiscible2ph")

#Set gravity
simulation.set_gravity(0)

#simulation.set_rock_volumetric_heat_capacity(rho_rock * cp_rock)

#Define the wells

# Initialize the regionalized values and distribute the domain
simulation.init(
    mesh=grid,
    #wells=make_wells,
    set_dirichlet_nodes=simulation.vertical_boundaries(grid),
    #set_global_rocktype=select_global_rocktype,
    cell_porosity=porosity(define_rock_mass_boundaries()[0],define_rock_mass_boundaries()[1]),
    cell_permeability=permeability(define_rock_mass_boundaries()[0],define_rock_mass_boundaries()[1]), 
    cell_thermal_conductivity=thconductivity(define_rock_mass_boundaries()[0],define_rock_mass_boundaries()[1])
)

# Setting up initial values
X0 = simulation.build_state(simulation.Context.liquid, p=pres, T=Tres)
simulation.all_states().set(X0)
simulation.dirichlet_node_states().set(X0)

# Define time steps parameters
simulation.standard_loop(
    initial_timestep=100 * day,
    final_time=30 * year,
    output_period=year,
)

simulation.postprocess()

I want to define two different regions with different phyiscal properties such as porosity, permeability,..... to run a numerical simulation in ComPASS software, however, the following error is coming and I don't know how can I sort it out. The error is Mesh vertices are not allocated." which I'm confused how to solve it.