You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Hello, first of all, congratulations for this amazing library @lululxvi !!
Well, I'm trying to replicate the methodology used in the article “https://doi.org/10.1016/j.jcp.2018.10.045” for the inverse problem for Navier Stokes equations. I'm using a file in CGNS format from CFD simulations. I'm trying to automate the code to recognize different boundary conditions, geometries, etc.
I have 2 questions:
1st Question
When I extract the variables and select points randomly from the total data set, it works, but when I try to separate the inlet, outlet, walls and interior, I can't get good convergence.
This is the code that gets good convergence and good results when I compare it to CFD files with the PINN result:
importh5pyimportnumpyasnp# Função para normalizar os dadosdefnormalize(data, min_val, max_val):
return (data-min_val) / (max_val-min_val)
# Função para desnormalizar os dadosdefdenormalize(data, min_val, max_val):
returndata* (max_val-min_val) +min_val# Função para extrair e normalizar todos os dados do arquivo CGNSdefextract_data(file_path):
withh5py.File(file_path, 'r') asf:
grid_coords=f['Base/Zone/GridCoordinates']
x_coords=grid_coords['CoordinateX/ data'][:]
y_coords=grid_coords['CoordinateY/ data'][:]
# Acessar soluções globais (velocidade e pressão)solution_data=f['Base/Zone/FlowSolution.N:1']
velocity_x=solution_data['VelocityX/ data'][:]
velocity_y=solution_data['VelocityY/ data'][:]
pressure=solution_data['Pressure/ data'][:]
# Armazenar os valores mínimos e máximos para desnormalizaçãomin_max_vals= {
'x_min': np.min(x_coords), 'x_max': np.max(x_coords),
'y_min': np.min(y_coords), 'y_max': np.max(y_coords),
'u_min': np.min(velocity_x), 'u_max': np.max(velocity_x),
'v_min': np.min(velocity_y), 'v_max': np.max(velocity_y),
'p_min': np.min(pressure), 'p_max': np.max(pressure),
}
# Normalizar os dadosx_coords_norm=normalize(x_coords, min_max_vals['x_min'], min_max_vals['x_max'])
y_coords_norm=normalize(y_coords, min_max_vals['y_min'], min_max_vals['y_max'])
velocity_x_norm=normalize(velocity_x, min_max_vals['u_min'], min_max_vals['u_max'])
velocity_y_norm=normalize(velocity_y, min_max_vals['v_min'], min_max_vals['v_max'])
pressure_norm=normalize(pressure, min_max_vals['p_min'], min_max_vals['p_max'])
returnx_coords_norm, y_coords_norm, velocity_x_norm, velocity_y_norm, pressure_norm, min_max_vals# Função para carregar dados de treinamento a partir do arquivo CGNSdefload_training_data_from_cgns(num):
# Carregar os dados diretamente do arquivo CGNSx_coords, y_coords, velocity_x, velocity_y, pressure, min_max_vals=extract_data('CGNSall.cgns')
N=x_coords.shape[0] # Número total de pontos disponíveis# Flatten para garantir que os dados estejam em formato (N, 1)x=x_coords.flatten()[:, None] # N x 1y=y_coords.flatten()[:, None] # N x 1u=velocity_x.flatten()[:, None] # N x 1v=velocity_y.flatten()[:, None] # N x 1p=pressure.flatten()[:, None] # N x 1# Concatenar todos os dados em um array para facilidade de manipulaçãodata=np.concatenate([x, y, u, v, p], axis=1) # N x 5# Escolher aleatoriamente `num` pontos de dadosidx=np.random.choice(N, num, replace=False)
selected_data=data[idx, :]
# Separar os dados selecionados em suas variáveisx_train=selected_data[:, 0:1] # Coordenada xy_train=selected_data[:, 1:2] # Coordenada yu_train=selected_data[:, 2:3] # Componente u da velocidadev_train=selected_data[:, 3:4] # Componente v da velocidadep_train=selected_data[:, 4:5] # Pressãoreturn [x_train, y_train, u_train, v_train, p_train, min_max_vals]
# Exemplo de usox_train, y_train, u_train, v_train, p_train, min_max_vals=load_training_data_from_cgns(num=1000)
# Caso seja necessário desnormalizar, podemos usar a função de desnormalizaçãox_desnorm=denormalize(x_train, min_max_vals['x_min'], min_max_vals['x_max'])
u_desnorm=denormalize(u_train, min_max_vals['u_min'], min_max_vals['u_max'])
importdeepxdeasddeimportnumpyasnp# Definir os parâmetros a serem identificadosC1=dde.Variable(0.0)
C2=dde.Variable(0.0)
# Define as equações de Navier-Stokes (para PDEs estacionárias)defNavier_Stokes_Equation(x, y):
u=y[:, 0:1]
v=y[:, 1:2]
p=y[:, 2:3]
# Derivadas de primeira ordemdu_x=dde.grad.jacobian(y, x, i=0, j=0)
du_y=dde.grad.jacobian(y, x, i=0, j=1)
dv_x=dde.grad.jacobian(y, x, i=1, j=0)
dv_y=dde.grad.jacobian(y, x, i=1, j=1)
dp_x=dde.grad.jacobian(y, x, i=2, j=0)
dp_y=dde.grad.jacobian(y, x, i=2, j=1)
# Derivadas de segunda ordem (laplaciano)du_xx=dde.grad.hessian(y, x, component=0, i=0, j=0)
du_yy=dde.grad.hessian(y, x, component=0, i=1, j=1)
dv_xx=dde.grad.hessian(y, x, component=1, i=0, j=0)
dv_yy=dde.grad.hessian(y, x, component=1, i=1, j=1)
# Continuidade e equações de momentocontinuity=du_x+dv_yx_momentum=C1* (u*du_x+v*du_y) +dp_x-C2* (du_xx+du_yy)
y_momentum=C1* (u*dv_x+v*dv_y) +dp_y-C2* (dv_xx+dv_yy)
return [continuity, x_momentum, y_momentum]
# Definir a geometria do domíniogeom=dde.geometry.Rectangle([0, 0], [2, 0.1])
# Carregar dados normalizados do arquivo CGNS e valores mínimos/máximos
[ob_x, ob_y, ob_u, ob_v, ob_p, min_max_vals] =load_training_data_from_cgns(num=1000)
ob_xy=np.hstack((ob_x, ob_y))
# Definir as condições de contorno observadas para u e vobserve_u=dde.icbc.PointSetBC(ob_xy, ob_u, component=0)
observe_v=dde.icbc.PointSetBC(ob_xy, ob_v, component=1)
# Definir os dados para o modelo PINNdata=dde.data.PDE(
geom,
Navier_Stokes_Equation,
[observe_u, observe_v],
num_domain=2000, # Pontos internos (collocation points)num_boundary=2000, # Pontos de contorno#anchors=ob_xy # Inclui as coordenadas observadas
)
# Configurar a rede neural (PINN)layer_size= [2] + [50] *8+ [3] # 2 entradas (x, y) e 3 saídas (u, v, p)activation="tanh"initializer="Glorot uniform"net=dde.maps.FNN(layer_size, activation, initializer)
# Criar o modelo PINNmodel=dde.Model(data, net)
# Callbacks para armazenar os resultados de C1 e C2fnamevar="variables.dat"variable=dde.callbacks.VariableValue([C1, C2], period=100, filename=fnamevar)
# Compilar o modelo e treinar com early stoppingmodel.compile("adam", lr=1e-5, external_trainable_variables=[C1, C2])
loss_history, train_state=model.train(iterations=50000, callbacks=[variable], display_every=500, disregard_previous_best=True)
# Salvar e plotar os resultadosdde.saveplot(loss_history, train_state, issave=True, isplot=True)
This is the convergence that seems good to me for this code:
However, when I try to divide the zones, I don't get good convergence because some residuals from the boundary conditions simply stabilize too early. This is the code with bad convergence that I want to improve:
importh5pyimportnumpyasnp# Função para normalizar os dados globalmentedefnormalize_global(data, min_val=None, max_val=None):
ifmin_valisNoneormax_valisNone:
min_val=np.nanmin(data)
max_val=np.nanmax(data)
return (data-min_val) / (max_val-min_val), min_val, max_val# Função para normalizar entre variáveis diferentesdefnormalize_between_vars(*arrays):
all_data=np.concatenate(arrays)
min_val=np.nanmin(all_data)
max_val=np.nanmax(all_data)
normalized_arrays= [(arr-min_val) / (max_val-min_val) forarrinarrays]
returnnormalized_arrays, min_val, max_val# Função para extrair e normalizar dados do arquivo CGNSdefextract_and_normalize(file_path):
withh5py.File(file_path, 'r') asf:
grid_coords=f['Base/Zone/GridCoordinates']
x_coords=grid_coords['CoordinateX/ data'][:]
y_coords=grid_coords['CoordinateY/ data'][:]
# Acessar soluções globais (velocidade e pressão)solution_data=f['Base/Zone/FlowSolution.N:1']
velocity_x=solution_data['VelocityX/ data'][:]
velocity_y=solution_data['VelocityY/ data'][:]
pressure=solution_data['Pressure/ data'][:]
# Normalização global para cada variávelx_coords_norm, x_min, x_max=normalize_global(x_coords)
y_coords_norm, y_min, y_max=normalize_global(y_coords)
velocity_x_norm, vx_min, vx_max=normalize_global(velocity_x)
velocity_y_norm, vy_min, vy_max=normalize_global(velocity_y)
pressure_norm, p_min, p_max=normalize_global(pressure)
zones_of_interest= ['interior-surface_body', 'inlet', 'outlet', 'walls']
zone_data= {}
forzoneinzones_of_interest:
try:
element_connectivity=f[f'Base/Zone/{zone}/ElementConnectivity/ data'][:]
unique_node_ids=np.unique(element_connectivity)
# Filtrar as coordenadas e soluções associadas a esses IDs de nószone_x_coords=x_coords_norm[unique_node_ids-1]
zone_y_coords=y_coords_norm[unique_node_ids-1]
zone_velocity_x=velocity_x_norm[unique_node_ids-1]
zone_velocity_y=velocity_y_norm[unique_node_ids-1]
zone_pressure=pressure_norm[unique_node_ids-1]
zone_data[zone] = {
'x_coords': zone_x_coords,
'y_coords': zone_y_coords,
'velocity_x': zone_velocity_x,
'velocity_y': zone_velocity_y,
'pressure': zone_pressure
}
exceptKeyErrorase:
print(f"Erro ao acessar a zona {zone}: {e}")
# Retorna apenas o necessárioreturnzone_data, {
'x_range': (x_min, x_max),
'y_range': (y_min, y_max),
'vx_range': (vx_min, vx_max),
'vy_range': (vy_min, vy_max),
'p_range': (p_min, p_max)
}
defselect_points_per_zone(zone_data, num_points_per_zone):
selected_zone_data= {}
forzone, datainzone_data.items():
total_points=len(data['x_coords'])
# Imprimir o número total de pontos disponíveis na zonaprint(f"Zona: {zone} tem {total_points} pontos disponíveis.")
# Verificar se o número de pontos na zona é menor que o solicitadoiftotal_points<num_points_per_zone:
print(f"Solicitado: {num_points_per_zone}, mas apenas {total_points} disponíveis. Selecionando todos.")
num_points=total_pointselse:
num_points=num_points_per_zone# Selecionar `num_points` pontos aleatóriosindices=np.random.choice(total_points, num_points, replace=False)
selected_zone_data[zone] = {
'x_coords': data['x_coords'][indices],
'y_coords': data['y_coords'][indices],
'velocity_x': data['velocity_x'][indices],
'velocity_y': data['velocity_y'][indices],
'pressure': data['pressure'][indices],
}
# Imprimir o número de pontos selecionadosprint(f"Zona: {zone} - {num_points} pontos selecionados\n")
returnselected_zone_dataimportdeepxdeasddeimportnumpyasnp# Configurar para usar float64 em todas as operaçõesdde.config.set_default_float("float64")
C1=dde.Variable(0.0, dtype="float64")
C2=dde.Variable(0.0, dtype="float64")
# Define as equações de Navier-Stokes (para PDEs estacionárias)defNavier_Stokes_Equation(x, y):
u=y[:, 0:1]
v=y[:, 1:2]
p=y[:, 2:3]
# Derivadas de primeira ordemdu_x=dde.grad.jacobian(y, x, i=0, j=0)
du_y=dde.grad.jacobian(y, x, i=0, j=1)
dv_x=dde.grad.jacobian(y, x, i=1, j=0)
dv_y=dde.grad.jacobian(y, x, i=1, j=1)
dp_x=dde.grad.jacobian(y, x, i=2, j=0)
dp_y=dde.grad.jacobian(y, x, i=2, j=1)
# Derivadas de segunda ordem (laplaciano)du_xx=dde.grad.hessian(y, x, component=0, i=0, j=0)
du_yy=dde.grad.hessian(y, x, component=0, i=1, j=1)
dv_xx=dde.grad.hessian(y, x, component=1, i=0, j=0)
dv_yy=dde.grad.hessian(y, x, component=1, i=1, j=1)
# Continuidade e equações de momentocontinuity=du_x+dv_yx_momentum=C1* (u*du_x+v*du_y) +dp_x-C2* (du_xx+du_yy)
y_momentum=C1* (u*dv_x+v*dv_y) +dp_y-C2* (dv_xx+dv_yy)
return [continuity, x_momentum, y_momentum]
# Função para desnormalizar os dadosdefdenormalize(data, min_val, max_val):
returndata* (max_val-min_val) +min_val# Extração e normalização dos dados do CGNSzone_data, min_max_vals=extract_and_normalize('CGNSall.cgns')
# Selecionar um número de pontos diferente para cada zonaselected_inlet_data=select_points_per_zone({'inlet': zone_data['inlet']}, 28)
selected_outlet_data=select_points_per_zone({'outlet': zone_data['outlet']}, 27)
selected_walls_data=select_points_per_zone({'walls': zone_data['walls']}, 161)
selected_interior_data=select_points_per_zone({'interior-surface_body': zone_data['interior-surface_body']}, 1000)
# Obter os dados selecionados de cada zonainlet_data=selected_inlet_data['inlet']
outlet_data=selected_outlet_data['outlet']
walls_data=selected_walls_data['walls']
interior_data=selected_interior_data['interior-surface_body']
inlet_xy=np.hstack((inlet_data['x_coords'][:, None], inlet_data['y_coords'][:, None]))
outlet_xy=np.hstack((outlet_data['x_coords'][:, None], outlet_data['y_coords'][:, None]))
walls_xy=np.hstack((walls_data['x_coords'][:, None], walls_data['y_coords'][:, None]))
interior_xy=np.hstack((interior_data['x_coords'][:, None], interior_data['y_coords'][:, None]))
# Obter as velocidades e pressão para cada zonainlet_u, inlet_v=inlet_data['velocity_x'], inlet_data['velocity_y']
outlet_u, outlet_v=outlet_data['velocity_x'], outlet_data['velocity_y']
walls_u, walls_v=walls_data['velocity_x'], walls_data['velocity_y']
interior_u, interior_v, interior_p=interior_data['velocity_x'], interior_data['velocity_y'], interior_data['pressure']
# Definir a geometria do domíniogeom=dde.geometry.Rectangle([0, 0], [2, 0.1])
# Definir as condições de contorno observadas para u e v (incluindo todas as zonas)observe_inlet_u=dde.icbc.PointSetBC(inlet_xy, inlet_u, component=0)
observe_inlet_v=dde.icbc.PointSetBC(inlet_xy, inlet_v, component=1)
observe_outlet_u=dde.icbc.PointSetBC(outlet_xy, outlet_u, component=0)
observe_outlet_v=dde.icbc.PointSetBC(outlet_xy, outlet_v, component=1)
observe_walls_u=dde.icbc.PointSetBC(walls_xy, walls_u, component=0)
observe_walls_v=dde.icbc.PointSetBC(walls_xy, walls_v, component=1)
observe_interior_u=dde.icbc.PointSetBC(interior_xy, interior_u, component=0)
observe_interior_v=dde.icbc.PointSetBC(interior_xy, interior_v, component=1)
# Definir os dados para o modelo PINN (incluindo as observações de todas as zonas)data=dde.data.PDE(
geom,
Navier_Stokes_Equation,
[observe_inlet_u, observe_inlet_v, observe_outlet_u, observe_outlet_v, observe_interior_u, observe_interior_v], #, observe_walls_u, observe_walls_v, num_domain=2000, # Pontos internos (collocation points)num_boundary=2000,
#anchors=interior_xy # Inclui as coordenadas observadas
)
# Configurar a rede neural (PINN)layer_size= [2] + [50] *8+ [3] # 2 entradas (x, y) e 3 saídas (u, v, p)activation="tanh"initializer="Glorot uniform"net=dde.maps.FNN(layer_size, activation, initializer)
# Criar o modelo PINNmodel=dde.Model(data, net)
# Callbacks para armazenar os resultados de C1 e C2 e early stoppingfnamevar="variables.dat"variable=dde.callbacks.VariableValue([C1, C2], period=100, filename=fnamevar)
early_stopping=dde.callbacks.EarlyStopping(min_delta=1e-4, patience=5000)
# Compilar o modelo e treinar com early stoppingmodel.compile("adam", lr=1e-5, external_trainable_variables=[C1, C2])
loss_history, train_state=model.train(iterations=500000, callbacks=[variable], display_every=1000, disregard_previous_best=True) #, early_stopping# Salvar e plotar os resultadosdde.saveplot(loss_history, train_state, issave=True, isplot=True)
This the convergence when I tried to separate al data in zones:
2nd Question
I don't care about the pressure results at the moment, but I would like to insert them in the future to get as good a resolution of the fluid for both velocities and pressures.
When it's working with the different zones, is it enough to teach with data in the bcs for pressures? I found somewhere that the convergence may not be the best and you have to use adam+L-FBGS.
I'd like to have a comment on this, because it is more data to the model get confuse.
Thank you very much for this incredible job with this library. Very useful!
reacted with thumbs up emoji reacted with thumbs down emoji reacted with laugh emoji reacted with hooray emoji reacted with confused emoji reacted with heart emoji reacted with rocket emoji reacted with eyes emoji
-
Hello, first of all, congratulations for this amazing library @lululxvi !!
Well, I'm trying to replicate the methodology used in the article “https://doi.org/10.1016/j.jcp.2018.10.045” for the inverse problem for Navier Stokes equations. I'm using a file in CGNS format from CFD simulations. I'm trying to automate the code to recognize different boundary conditions, geometries, etc.
I have 2 questions:
1st Question
When I extract the variables and select points randomly from the total data set, it works, but when I try to separate the inlet, outlet, walls and interior, I can't get good convergence.
Data file : CGNSall.zip
This is the code that gets good convergence and good results when I compare it to CFD files with the PINN result:
This is the convergence that seems good to me for this code:
However, when I try to divide the zones, I don't get good convergence because some residuals from the boundary conditions simply stabilize too early. This is the code with bad convergence that I want to improve:
This the convergence when I tried to separate al data in zones:
2nd Question
I don't care about the pressure results at the moment, but I would like to insert them in the future to get as good a resolution of the fluid for both velocities and pressures.
When it's working with the different zones, is it enough to teach with data in the bcs for pressures? I found somewhere that the convergence may not be the best and you have to use adam+L-FBGS.
I'd like to have a comment on this, because it is more data to the model get confuse.
Thank you very much for this incredible job with this library. Very useful!
Beta Was this translation helpful? Give feedback.
All reactions