Comunicación Fortran - C++

Para realizar una simulación realista del efecto que provoca el control sobre el modelo físico es necesario aplicar la actuación en simultaneo con la resolución del sistema. Para ello se plantea un esquema de comunicación entre el Code Saturne y nuestro programa en C++ donde podamos definir el valor de actuación para cada paso de la simulación.
Se cuenta con varias posibilidades concretas de comunicación IPC: sockets, pipes, archivos. En este caso se opta por la comunicación mediante archivos FIFO (es decir, implementando el concepto de pipes mediante archivos virtuales) donde un programa lee del archivo mientras el otro escribe en el mismo.
Para analizar la viabilidad se decide realizar un programa simple donde dos procesos se comuniquen mutuamente. En este caso contamos con un proceso que lee e imprime los valores de actuación a utilizar (hecho en Fortran 90) y otro proceso que genera y escribe los mismos (C++).

Lector (Fortran 90) - lector.f90

program readcsv
  implicit none
  call readAndPrint()
end program readcsv

subroutine readAndPrint()
  integer status
  double precision value

  open(unit=99, file="act1.csv", action="read")
  do
    read(99, *, IOSTAT=status) value
    if (status > 0) then
      !Error
      write (102,*) 'An error has been detected in the file'
      exit
    else if (status < 0) then
      !EOF
      exit
    else
      print *, 'Value:',value
    end if
  end do
  close(99)
end subroutine

Escritor (C++) - escritor.cpp

#include <iostream>
#include <fstream>
#include <ctime>

int main(int argc, char* argv[]) {
        std::ofstream arch("act1.csv", std::fstream::app);
        for (int i = 0; i < 2; ++i) {
                arch << 1.0 + i / 100.0 << std::endl;
                sleep(2);
        }
        return 0;
} 

Para ejecutar el programa es necesario abrir un archivo FIFO del sistema operativo. Esto es muy sencillo en sistemas Linux/Unix:

mkfifo act1.csv

La secuencia real de compilación y ejecución es la siguiente:

$ gfortran -o lector lector.f90
$ g++ -o escritor escritor.cpp
$ mkfifo act1.csv
$ ./lector &
[1] 892 

$ ./escritor

 Value:   1.0000000000000000     
 Value:   1.0100000000000000     
[1]+  Done                    ./lector

Reducción de Variables y agregado de Trazadores

Sensores

Luego de las primeras pruebas con el algoritmo ARX propio del laboratorio, y teniendo en cuenta los pobres resultados de performance y precisión para grandes cantidades de variables, se decide simplificar el problema de estudio a un conjunto reducido de sensores.
Como primer aproximación se plantea la utilización de 4 únicos puntos sensores y la medición de la presión en esos puntos.

Los puntos se encuentran simétricos al eje X por lo que se espera una función espejada en cada par de los mismos. Los puntos colocados en la segunda línea de sensores deberían mantener cierta similitud con los de la primer línea pero con un decaimiento en la intensidad debido a la difusión más un retardo en el tiempo por la velocidad del escurrimiento.

Trazadores

Si bien esta información se puede obtener en el simulador CodeSaturne, es necesario idear un método de procesamiento que permita calcularla de las imágenes a obtener en el experimento real. Por tal motivo se decide agregar trazadores a la simulación que serán implementados con fuentes de humo (u otras sustancias que queden en suspención sobre el fluido) al montar el experimento.
A tal fin se agregan escalares a la simulación del Code Saturne como se observa a continuación:

ARX - Resultados del Algoritmo con Información Real

Luego de escribir el algoritmo ARX que soporte grandes vólumenes de información se decide ponerlo en práctica con los datos experimentales obtenidos de la simulación del Code Saturne.
De esta forma se pretende obtener resultados más cercanos a los esperados de la experiencia real y comprobar que los algoritmos escritos hasta el momento sean efectivos y eficaces.
Se practican varias pruebas con los valores de vorticidad obtenidos de la simulación obtenidos mediante una grilla definida en la zona de la estela del cilindro. Se colocan varias líneas de obtención de información y por cada línea se relevan aprox. 50 puntos en 500 pasos de tiempo a lo largo de 20 líneas.

Los resultados obtenidos con el total de la información no fueron satisfactorios. Las simulaciones divergen y se demoran mucho tiempo.
A fin de determinar un buen número de muestras a utilizar se realizan distintas comparaciones planteando configuraciones de: cantidad de líneas utilizadas, cantidad de muestras por línea, cantidad de pasos de simulación a predecir (antes de poder relevar los valores reales y poder realimentar al predictor). Los resultados fueron positivos siempre que se utilizaran menos de 200 variables de entrada. Se puede ver el porcentaje de ajuste para las distintas configuraciones en el siguiente gráfico:

ARX - Algoritmo propio para identificación del sistema

El algoritmo provisto por Matlab demostró funcionar a la perfección con las series de ejemplo pero al intentar utilizarlo con volúmenes de datos más grandes los tiempos de cómputo crecen exponencialmente. Como agregado, el sistema final debe programarse en C/C++/CUDA donde la utilización del toolkit de Matlab es poco práctica.
En consecuencia se decide implementar la identificación del modelo ARX. A tal fin se analizan ciertas publicaciones de Lennart Ljung (autor del toolkit de Matlab) y se consigue un algoritmo que entrega resultados similares en menor tiempo.
Para comprobar la bondad del algoritmo creado se realizan comparaciones contra los modelos obtenidos con el toolkit de Matlab para las series 'Motorized Camera' y 'Missile' previamente estudiadas.
Se utiliza la función compare del toolkit para verificar el ajuste con muy buenos resultados. Las siguientes imágenes muestran la comparación gráfica entre los resultados obtenidos por Matlab (llamado model) y por el algoritmo basado en los escritos de Ljung (llamado modelLjung):

Serie 'Motorized Camera'. Comparación entre algoritmos de identificación.

Serie 'Missile'. Comparación entre algoritmos de identificación.

ARX - Identificación del Sistema

El modelo físico a estudiar debe estar representado en el software de control de forma tal que sea posible predecir los futuros estados y determinar cual es la mejor actuación para poder conseguir el resultado deseado.
Existen varias maneras de modelar numéricamente una experiencia física. En este caso optamos por un modelo de caja negra llamado Auto Regressive with eXogenous input (ARX). En este modelo, es necesario definir grandes matrices de coeficientes que permitan, dado un vector de estado con información de pasos previos, obtener una predicción a futuro del sistema. Como todo modelo de caja cerrada, los coeficientes necesarios son definidos 'por inspección' o con métodos de identificación en base a mediciones previas.

La ecuación que define el modelo ARX más simple, una entrada y una salida (SISO), es la siguiente:

donde:

y: variable de salida del sistema
u: variable de entrada del sistema
a y b: coeficientes asociados

En caso de tener un sistema con entradas y salidas múltiples (MIMO), las variables u e y se transforman en vectores (un componente de u por cada entrada, un componente de y por cada salida). Sin cambios en la ecuación, se la puede extender a la forma matricial para representar ese esquema.
Por último, los coeficientes a y b se pueden agrupar utilizando una notación reducida:

Donde se definen las matrices de polinomios A y B según:

Siendo Z el operador de corrimiento:

En todos los casos se debe entender la existencia de un término de error en cualquiera de los lados de la igualdad que corrige las incertezas que provoca el modelo a caja cerrada. Con esto en mente, se puede entender el modelo con el siguiente esquema:

Matlab posee un toolkit de identificación de sistemas que entrega el modelo ARX computado a partir de una serie de muestras. Para probar el toolkit utilizamos series predefinidas en Matlab llamadas 'Motorized Camera' y 'Missile' que incluyen datos reales de varios sensores.
Se puede obtener más información sobre el algoritmo y las funciones utilizadas desde aquí.

% REAL DATA 1: MISSILE
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'missiledata.mat'));
% o bien: 

% load motorizedcamera
data = iddata(y, u);

p = 2;

outputs = size(y,2);
inputs = size(u,2);
model = arx(data, [p*ones(outputs, outputs), p*ones(outputs, inputs), ones(outputs, inputs)]);

Video de la Simulación con Actuadores Programados

Gracias al código Fortran agregado a los loops de simulación del CodeSaturne, pudimos agregar actuadores de plasma a los resultados anteriores.
De esta forma simulamos la actuación en dos sectores del cilindro, realizando pulsaciones de plasma con ángulo de salida tangente para ver la respuesta del sistema. Se espera que las turbulencias del cilindro desaparezcan con el plasma aplicado y vuelvan a generarse tan pronto se retire el efecto del actuador.
Durante esta simulación se agregaron algunos esquemas de actuación periódicos, oscilantes y asiméticos para poder relevar información que será usada en la identificación del sistema mediante el modelo ARX.

Usando variables 'static' en Fortran para cargar un CSV sólo una vez

Para completar la idea de utilizar un vector desde un archivo dentro del código de code saturne es importante dar ciertas garantías de performance. La subrutina a utilizar se encuentra dentro del archivo usclim.f90 que es utilizado de forma repetitiva durante la simulación.

Una forma de garantizar que la carga del CSV se realiza una sóla vez es utilizando código global y variables globales. Lamentablemente, el archivo disponible no cuenta con las secciones comunes con lo que se supone una composición de ese código dentro del programa de simulación general.

Otra alternativa es utilizar variables en el data segment persistentes entre llamadas a la subrutina. Optamos por ese camino y utilizamos la palabra reservada SAVE para realizar un breve programa de prueba:

program read_csv_once
  implicit none
 
  call readAndPrint(1000)
  call readAndPrint(100)

end program read_csv_once

subroutine readAndPrint(test_offset)
  integer:: test_offset
  integer, parameter:: values_qty = 2
  real, save, pointer:: pValues(:)
  logical, save:: pValues_initialized= .false.

  if (.not.pValues_initialized) then
    allocate(pValues(values_qty))
    open(unit=99, file="act1.csv", action="read")
    read(99, *) pValues
    close(99)
    pValues = pValues + test_offset
    pValues_initialized = .true.
  endif

  write(*,*) "Values from file:", pValues
  write(*,*) "Values Quantity:", values_qty
end subroutine

Nuevamente tenemos una cantidad fija de valores a leer del CSV dado que por el momento es un número conocido.
Notar que el puntero estático pValues NO posee un valor por defecto definido, es por ese motivo que se utiliza pValues_initialized para controlar si se hizo la lectura del CSV o no.
El resultado de correr este programa de pruebas es:

 Values from file:   1003.6295       1003.6688    
 Values Quantity:           2
 Values from file:   1003.6295       1003.6688    
 Values Quantity:           2

Demostrando que la inicialización se realizó una vez.