ComPASS模拟定义双区域物性时遇“Mesh vertices未分配”错误求助

问题描述 #

尝试在ComPASS中定义两个具有不同孔隙度、渗透率等物理属性的区域进行数值模拟时，运行代码出现错误：Mesh vertices are not allocated。

错误代码 #

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()

错误原因分析 #

1. 网格初始化时机错误：simulation.compute_global_cell_centers()需要在simulation.init()之后调用，只有初始化后网格顶点和单元信息才会被分配，提前调用会触发Mesh vertices are not allocated错误。
2. 变量作用域问题：porosity、permeability等函数中直接使用cell_centers变量，但该变量仅在define_rock_mass_boundaries()内部定义，会引发NameError。
3. 重复计算浪费资源：多次调用define_rock_mass_boundaries()会重复计算单元中心和区域索引，效率低下且可能引发意外问题。
4. 物理属性数组形状错误：np.ones(cell_centers.shape)创建的是三维数组（每个单元3个坐标），但ComPASS要求物理属性是一维数组（每个单元对应一个值）。
5. 区域坐标范围错误：网格原点Ox=-1500，Lx=4000，实际x坐标范围是-1500到2500，原代码中0到4000的范围会导致大部分单元无法被匹配。

修复后的代码 #

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

# 设置输出信息
ComPASS.set_output_directory_and_logfile(__file__)

# 参数定义
pres = 30. * MPa            # 初始储层压力
Tres = degC2K(200.)         # 初始储层温度（转换为开尔文）
K_reservoir = 2             # 体积热导率 W/m/K
g = 9.81                    # 重力加速度 m/s^2

# 网格参数
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)

# 加载物理模型
simulation = ComPASS.load_physics("immiscible2ph")
simulation.set_gravity(0)

# 初始化模拟（先完成网格初始化）
simulation.init(
    mesh=grid,
    set_dirichlet_nodes=simulation.vertical_boundaries(grid),
)

# 计算单元中心并划分区域
cell_centers = simulation.compute_global_cell_centers()
# 基于实际网格x范围（-1500到2500）划分区域，示例以x=500为界
max_x_cox = 500
min_x_cct = 500
max_x_cct = Ox + Lx

cox_indices = np.where((cell_centers[:, 0] >= Ox) & (cell_centers[:, 0] < max_x_cox))[0]
cct_indices = np.where((cell_centers[:, 0] >= min_x_cct) & (cell_centers[:, 0] <= max_x_cct))[0]

print("cell_centers shape:", cell_centers.shape)
print("cct indices count:", len(cct_indices))
print("cox indices count:", len(cox_indices))

# 定义物理属性函数
def get_porosity():
    porosity_grid = np.ones(nx*ny*nz)
    porosity_grid[cox_indices] = 0.15
    porosity_grid[cct_indices] = 0.3
    return porosity_grid

def get_permeability():
    permeability_grid = np.ones(nx*ny*nz)
    permeability_grid[cox_indices] = 1e-14
    permeability_grid[cct_indices] = 1e-18
    return permeability_grid

def get_thermal_conductivity():
    thcond_grid = np.ones(nx*ny*nz)
    thcond_grid[cox_indices] = 2.1
    thcond_grid[cct_indices] = 4
    return thcond_grid

# 设置单元物理属性
simulation.set_cell_porosity(get_porosity())
simulation.set_cell_permeability(get_permeability())
simulation.set_cell_thermal_conductivity(get_thermal_conductivity())

# 设置初始状态
X0 = simulation.build_state(simulation.Context.liquid, p=pres, T=Tres)
simulation.all_states().set(X0)
simulation.dirichlet_node_states().set(X0)

# 运行模拟循环
simulation.standard_loop(
    initial_timestep=100 * day,
    final_time=30 * year,
    output_period=year,
)

simulation.postprocess()

关键修改说明 #

• 调整初始化顺序：先调用simulation.init()完成网格初始化，再计算单元中心和划分区域，避免“Mesh vertices are not allocated”错误。
• 修正变量作用域：将cell_centers提升到全局作用域，确保属性函数可以访问；或者直接在函数内部使用已计算的索引。
• 统一区域计算：只计算一次区域索引，避免重复调用带来的性能损耗和错误。
• 修正数组形状：物理属性数组使用np.ones(nx*ny*nz)创建一维数组，符合ComPASS的要求。
• 修正区域坐标范围：基于网格实际的x坐标范围划分区域，确保单元能被正确匹配。

内容的提问来源于stack exchange，提问作者Afshin Davarpanah