text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
/* inet_aton -> per utilizzarla in gcc: __USE_MISC */
#define __USE_MISC
#include <sys/types.h>
#include <sys/socket.h>
#include <netinet/in.h>
#include <arpa/inet.h>
#include <pthread.h>
#include "buffer.h"
#include "list.h"
#include "router.h"
#include <gsl/gsl_matrix.h>

#define DIM_BUFFER 1000
#define POISON_ID -1

struct msg_t {
    int msg_id;
    int msg_time;
    float msg_x;
    float msg_z;
};

struct routing_table_t {
    int dimension;
    list list_table;
};

struct router_t {
    int router_id;
    struct routing_table_t routing_table;
    pthread_mutex_t router_mutex;
};

struct parametri_accepter_t {
    int accepter_agents_number;
    char *accepter_address;
    int accepter_port;
    list accepter_lista_buffer;
    int accepter_in_out;
};

struct parametri_dispatcher_t {
    buffer dispatcher_buffer_in;
    list dispatcher_lista_buffer_out;
};

struct parametri_receiver_t {
    int receiver_id;
    int receiver_socket;
    buffer receiver_buffer;
};

struct parametri_sender_t {
    int sender_id;
    int sender_socket;
    buffer sender_buffer;
};

typedef struct msg_t *msg;

typedef struct parametri_accepter_t *parametri_accepter;
typedef struct parametri_dispatcher_t *parametri_dispatcher;
typedef struct parametri_receiver_t *parametri_receiver;
typedef struct parametri_sender_t *parametri_sender;

static pthread_mutex_t router_class_mutex;
static int routerObjectsNumber = 0;

static int numero_receiver_vivi = 0;
static int numero_dispatcher_vivi = 0;
static int numero_sender_vivi = 0;

static void initClassRouter() {
    pthread_mutex_init(&(router_class_mutex), NULL);
}

static void cleanClassRouter() {
    pthread_mutex_destroy(&(router_class_mutex));

}

router allocRouter() {

    router router_M_;

    if (routerObjectsNumber == 0) initClassRouter();

    routerObjectsNumber++;

    router_M_ = (router) malloc(sizeof (struct router_t));

    pthread_mutex_init(&(router_M_->router_mutex), NULL);

    return router_M_;

}

void freeRouter(router _F_router) {

    pthread_mutex_destroy(&(_F_router->router_mutex));

    free(_F_router);

    _F_router = NULL;

    routerObjectsNumber--;

    if (routerObjectsNumber == 0) cleanClassRouter();
}

static int listenFrom(char * _ip_address_local, int _ip_port_local) {

    int socket_M_;

    struct sockaddr_in sockaddr_in_local_;

    memset(&sockaddr_in_local_, 0, sizeof (sockaddr_in_local_));
    sockaddr_in_local_.sin_family = AF_INET;
    sockaddr_in_local_.sin_port = htons(_ip_port_local);
    inet_aton(_ip_address_local, &(sockaddr_in_local_.sin_addr));

    if ((socket_M_ = socket(AF_INET, SOCK_STREAM, 0)) == -1) {
        perror("listenFrom: socket() Error\n");
        exit(1);
    }

    if (bind(socket_M_, (struct sockaddr *) & sockaddr_in_local_, sizeof (struct sockaddr_in)) == -1) {
        perror("listenFrom: bind() Error\n");
        exit(1);
    }

    if (listen(socket_M_, 1) == -1) {
        perror("listenFrom: listen() Error\n");
        exit(1);
    }

    return socket_M_;

}

static void sendMsg(int _socket, msg _msg) {

    if (send(_socket, _msg, sizeof (struct msg_t), 0) == -1) {
        perror("sendMsg: send() Error\n");
        exit(1);
    };
    /*
        int net_id;
        int net_x;
        int net_z;

        net_id = htonl(_msg->msg_id);
        if (send(_socket, &net_id, sizeof (net_id), 0) == -1) {
            perror("sendMsg: send() Error\n");
            exit(1);
        };

        net_x = htonl(_msg->msg_x);
        if (send(_socket, &net_x, sizeof (net_x), 0) == -1) {
            perror("sendMsg: send() Error\n");
            exit(1);
        };

        net_z = htonl(_msg->msg_z);
        if (send(_socket, &net_z, sizeof (net_z), 0) == -1) {
            perror("sendMsg: send() Error\n");
            exit(1);
        };
     */
}

static msg recvMsg(int _socket) {

    msg _M_msg = (struct msg_t *) malloc(sizeof (struct msg_t));

    if (recv(_socket, _M_msg, sizeof (struct msg_t), 0) == -1) {
        perror("recvMsg: recv() Error\n");
        exit(1);
    }
    /*
        int net_id;
        int net_x;
        int net_z;

        msg _M_msg = (struct msg_t *) malloc(sizeof (struct msg_t));

        if (recv(_socket, &net_id, sizeof (net_id), 0) == -1) {
            perror("recvMsg: recv() Error\n");
            exit(1);
        }
        _M_msg->msg_id = ntohl(net_id);

        if (recv(_socket, &net_x, sizeof (net_x), 0) == -1) {
            perror("recvMsg: recv() Error\n");
            exit(1);
        }
        _M_msg->msg_x = ntohl(net_x);

        if (recv(_socket, &net_z, sizeof (net_z), 0) == -1) {
            perror("recvMsg: recv() Error\n");
            exit(1);
        }
        _M_msg->msg_z = ntohl(net_z);
     */
    return _M_msg;
}

static parametri_receiver allocParamReceiver(int _receiver_socket, buffer _receiver_buffer) {

    parametri_receiver _M_parametri_receiver = (parametri_receiver) malloc(sizeof (struct parametri_receiver_t));

    _M_parametri_receiver->receiver_socket = _receiver_socket;
    _M_parametri_receiver->receiver_buffer = _receiver_buffer;

    return _M_parametri_receiver;
}

/*
 * Quando il Receiver riceve dall'agent un messaggio con id = POISON_ID allora:
 * - lo inoltra nel buffer
 * - chiude la connessione con l'agent
 * - fa il fflush del file di log e lo chiude
 * - decrementa numero_receiver_vivi
 * - termina
 */

static void *runReceiver(void *_parametri_receiver) {

    parametri_receiver _parametri = (parametri_receiver) _parametri_receiver;

    int _receiver_socket = _parametri->receiver_socket;
    buffer _receiver_buffer = _parametri->receiver_buffer;

    msg _msg;
    int _receiver_id;

    FILE *receiver_log;

    char receiver_log_file_name[FILENAME_MAX];

    /*
     * Riceve il messaggio di presentazione dall'agent
     * - imposta l'id del receiver
     * - imposta l'id del _receiver_buffer
     * - reinvia all'agent il messaggio ricevuto
     */

    _msg = recvMsg(_receiver_socket);
    _receiver_id = _msg->msg_id;
    setBufferID(_receiver_buffer, _receiver_id);

    /*
     * Apre il file di log
     */

    sprintf(receiver_log_file_name, "receiver_%d.log", _receiver_id);

    receiver_log = fopen(receiver_log_file_name, "w");

    if (receiver_log == NULL) {
        perror("Errore apertura file: \"receiver.log\"\n");
        exit(1);
    };

    setbuf(receiver_log, NULL);

    fprintf(receiver_log, "R%d: Open log\n", _receiver_id);

    /*sleep(1);*/

    fprintf(receiver_log, "R%d: Parametri invocazione:\n", _receiver_id);

    fprintf(receiver_log, "R%d: Wait: rcv\n", _receiver_id);

    fprintf(receiver_log, "R%d: > (%d,%d,%f,%f)\n", _receiver_id, _msg->msg_id, _msg->msg_time, _msg->msg_x, _msg->msg_z);

    fprintf(receiver_log, "R%d: Wait: rcv\n", _receiver_id);

    /*
     * Riceve un messaggio e lo pone nel buffer finchè non riceve il POISON_ID
     */

    _msg = recvMsg(_receiver_socket);

    fprintf(receiver_log, "R%d: > (%d,%d,%f,%f)\n", _receiver_id, _msg->msg_id, _msg->msg_time, _msg->msg_x, _msg->msg_z);

    fprintf(receiver_log, "R%d: Loop:\n", _receiver_id);

    fprintf(receiver_log, "R%d: > while(msg != POISON_ID)\n", _receiver_id);

    fprintf(receiver_log, "R%d: --- Loop Start ---\n", _receiver_id);

    while (_msg->msg_id != POISON_ID) {

        putInBuffer_blocking(_receiver_buffer, _msg);

        fprintf(receiver_log, "R%d: >>> putInBuffer(%d,%d,%f,%f)\n", _receiver_id, _msg->msg_id, _msg->msg_time, _msg->msg_x, _msg->msg_z);

        fprintf(receiver_log, "R%d: >>> Wait: rcv\n", _receiver_id);

        _msg = recvMsg(_receiver_socket);
    };

    fprintf(receiver_log, "R%d: ---- Loop End ----\n", _receiver_id);

    fprintf(receiver_log, "R%d: Signal:\n", _receiver_id);

    putInBuffer_blocking(_receiver_buffer, _msg);

    fprintf(receiver_log, "R%d: > putInBuffer(%d,%d,%f,%f) \n", _receiver_id, _msg->msg_id, _msg->msg_time, _msg->msg_x, _msg->msg_z);

    fprintf(receiver_log, "R%d: Socket:\n", _receiver_id);

    /*
     * Chiude la connessione con l'agent
     */

    send(_receiver_socket, "", 1, 0);

    close(_receiver_socket);

    fprintf(receiver_log, "R%d: > close receiver socket\n", _receiver_id);

    fprintf(receiver_log, "R%d: Memoria\n", _receiver_id);

    free(_parametri_receiver);

    fprintf(receiver_log, "R%d: > parametri\n", _receiver_id);

    /*
     * fflush e chiusura del logfile
     */

    fprintf(receiver_log, "R%d: Close log\n", _receiver_id);

    fflush(receiver_log);

    fclose(receiver_log);

    /* Decrementa numero_receiver_vivi e termina*/

    numero_receiver_vivi--;

    return NULL;
}

static parametri_sender allocParamSender(int _sender_socket, buffer _sender_buffer) {

    parametri_sender _M_parametri_sender = (parametri_sender) malloc(sizeof (struct parametri_sender_t));

    _M_parametri_sender->sender_socket = _sender_socket;
    _M_parametri_sender->sender_buffer = _sender_buffer;

    return _M_parametri_sender;
}

/*
 * Quando il Sender prende dal buffer un messaggio con id = POISON_ID
 * vuol dire che quello è l'ultimo messaggio e il buffer non ne contiene altri:
 * - attende che i dispatcher siano morti
 * - elimina il sender buffer associato
 * - chiude la connessione con l'agent
 * - fa il fflush del file di log e lo chiude
 * - decrementa numero_sender_buffer
 * - termina
 */

static void *runSender(void *_parametri_sender) {

    parametri_sender _parametri = (parametri_sender) _parametri_sender;

    int _sender_socket = _parametri->sender_socket;
    buffer _sender_buffer = _parametri->sender_buffer;

    int _sender_id;

    msg _msg;

    FILE *sender_log;

    char sender_log_file_name[FILENAME_MAX];

    /*
     * Riceve il messaggio di presentazione dall'agent:
     * - imposta l'id del sender
     * - imposta l'id del _sender_buffer
     * - reinvia all'agent il messaggio ricevuto
     */

    _msg = recvMsg(_sender_socket);
    _sender_id = _msg->msg_id;
    setBufferID(_sender_buffer, _sender_id);

    /*
     * Apre il file di log
     */

    sprintf(sender_log_file_name, "sender_%d.log", _sender_id);

    sender_log = fopen(sender_log_file_name, "w");

    if (sender_log == NULL) {
        perror("Errore apertura file: \"sender.log\"\n");
        exit(1);
    };

    setbuf(sender_log, NULL);

    fprintf(sender_log, "S%d: Open log\n", _sender_id);

    /*sleep(1);*/

    fprintf(sender_log, "S%d: Parametri di invocazione:\n", _sender_id);

    fprintf(sender_log, "S%d: Wait: rcv\n", _sender_id);

    fprintf(sender_log, "S%d: > (%d,%d,%f,%f)\n", _sender_id, _msg->msg_id, _msg->msg_time, _msg->msg_x, _msg->msg_z);

    fprintf(sender_log, "S%d: Wait: getFromBuffer\n", _sender_id);

    /*
     * Riceve un messaggio e lo pone nel buffer finchè non riceve il POISON_ID
     */

    _msg = getFromBuffer_blocking(_sender_buffer);

    fprintf(sender_log, "S%d: > (%d,%d,%f,%f)\n", _sender_id, _msg->msg_id, _msg->msg_time, _msg->msg_x, _msg->msg_z);

    fprintf(sender_log, "S%d: Loop:\n", _sender_id);

    fprintf(sender_log, "S%d: > while(msg != POISON_ID)\n", _sender_id);

    fprintf(sender_log, "S%d: --- Loop Start ---\n", _sender_id);

    while (_msg->msg_id != POISON_ID) {

        sendMsg(_sender_socket, _msg);

        fprintf(sender_log, "S%d: >>> snd (%d,%d,%f,%f)\n", _sender_id, _msg->msg_id, _msg->msg_time, _msg->msg_x, _msg->msg_z);

        fprintf(sender_log, "S%d: >>> Wait: getFromBuffer\n", _sender_id);

        _msg = getFromBuffer_blocking(_sender_buffer);
    };

    fprintf(sender_log, "S%d: ---- Loop End ----\n", _sender_id);

    fprintf(sender_log, "S%d: Wait:\n", _sender_id);

    /*
     * Attende che tutti i dispatcher siano morti
     */

    fprintf(sender_log, "S%d: > while(dispatcher > 0)\n", _sender_id);

    while (numero_dispatcher_vivi > 0);

    fprintf(sender_log, "S%d: Socket:\n", _sender_id);

    /*
     * Chiude la connessione con l'agent
     */

    send(_sender_socket, "", 1, 0);

    close(_sender_socket);

    fprintf(sender_log, "S%d: > close sender socket\n", _sender_id);

    fprintf(sender_log, "S%d: Memoria:\n", _sender_id);

    /*
     * Elimina il receiver buffer associato
     */

    freeBuffer(_sender_buffer);

    fprintf(sender_log, "S%d: > free(sender_buffer)\n", _sender_id);

    free(_parametri_sender);

    fprintf(sender_log, "S%d: > parametri\n", _sender_id);

    /*
     * fflush e chiusura del logfile
     */

    fprintf(sender_log, "S%d: Close log\n", _sender_id);

    fflush(sender_log);

    fclose(sender_log);

    /* Decrementa numero_sender_vivi e termina*/

    numero_sender_vivi--;

    return NULL;
}

static parametri_accepter allocParamAccepter(int _accepter_agents_number, char *_accepter_address, int _accepter_port, list _accepter_lista_buffer, int _accepter_in_out) {

    parametri_accepter _M_parametri_accepter = (parametri_accepter) malloc(sizeof (struct parametri_accepter_t));

    _M_parametri_accepter->accepter_agents_number = _accepter_agents_number;
    _M_parametri_accepter->accepter_address = _accepter_address;
    _M_parametri_accepter->accepter_port = _accepter_port;
    _M_parametri_accepter->accepter_lista_buffer = _accepter_lista_buffer;
    _M_parametri_accepter->accepter_in_out = _accepter_in_out;

    return _M_parametri_accepter;

}

/*
 * Quando l'accepter riceve _accepter_agents_number connessioni in input e in output
 * vuol dire che tutti gli agents si sono collegati:
 * - chiude le due porte in listening
 * - fa il fflush del file di log e lo chiude
 * - termina
 */

static void *runAccepter(void *_parametri_accepter) {

    parametri_accepter _parametri = (parametri_accepter) _parametri_accepter;

    int _accepter_agents_number;
    char *_accepter_address;
    int _accepter_port;
    list _accepter_lista_buffer;
    int _accepter_in_out;

    int listening_socket;
    int communication_socket;

    FILE *accepter_log;

    char *accepter_log_file_name;
    char *in_out;
    char *accepter_error;

    /* Apre il file di log */

    _accepter_in_out = _parametri->accepter_in_out;

    if (_accepter_in_out == 0) {
        accepter_log_file_name = "accepter_in.log";
        accepter_error = "Errore apertura file: \"accepter_in.log\"\n";
        in_out = "AI";
    } else {
        accepter_log_file_name = "accepter_out.log";
        accepter_error = "Errore apertura file: \"accepter_out.log\"\n";
        in_out = "AO";
    };

    accepter_log = fopen(accepter_log_file_name, "w");

    if (accepter_log == NULL) {
        perror(accepter_error);
        exit(1);
    };

    setbuf(accepter_log, NULL);

    fprintf(accepter_log, "%s: Open log\n", in_out);

    fprintf(accepter_log, "%s: Parametri di invocazione:\n", in_out);

    _accepter_agents_number = _parametri->accepter_agents_number;

    fprintf(accepter_log, "%s: > %d\n", in_out, _accepter_agents_number);

    _accepter_address = _parametri->accepter_address;

    fprintf(accepter_log, "%s: > %s\n", in_out, _accepter_address);

    _accepter_port = _parametri->accepter_port;

    fprintf(accepter_log, "%s: > %d\n", in_out, _accepter_port);

    _accepter_lista_buffer = _parametri->accepter_lista_buffer;

    fprintf(accepter_log, "%s: > %d\n", in_out, _accepter_in_out);

    fprintf(accepter_log, "%s: Socket:\n", in_out);

    fprintf(accepter_log, "%s: > listenFrom(%s,%d)\n", in_out, _accepter_address, _accepter_port);

    listening_socket = listenFrom(_accepter_address, _accepter_port);

    if (_accepter_in_out == 0) {

        fprintf(accepter_log, "%s: Loop:\n", in_out);

        fprintf(accepter_log, "%s: > while(receiver < agents)\n", in_out);

        fprintf(accepter_log, "%s: ---------- Loop Start ----------\n", in_out);

        while (numero_receiver_vivi < _accepter_agents_number) {

            buffer _M_buffer = allocBuffer(DIM_BUFFER);

            pthread_t receiver_;

            if ((communication_socket = accept(listening_socket, NULL, NULL)) == -1) {
                perror("Accepter accept() Error\n");
                exit(1);
            }

            fprintf(accepter_log, "%s: >>> accept connection\n", in_out);

            addElementToList(_accepter_lista_buffer, _M_buffer);

            numero_receiver_vivi++;

            pthread_create(&receiver_, NULL, &runReceiver, allocParamReceiver(communication_socket, _M_buffer));
        }
        fprintf(accepter_log, "%s: ----------- Loop End -----------\n", in_out);

    } else {

        fprintf(accepter_log, "%s: Loop:\n", in_out);

        fprintf(accepter_log, "%s: > while(sender < agents)\n", in_out);

        fprintf(accepter_log, "%s: --------- Loop Start ---------\n", in_out);

        while (numero_sender_vivi < _accepter_agents_number) {

            buffer _M_buffer = allocBuffer(DIM_BUFFER);

            pthread_t sender_;

            if ((communication_socket = accept(listening_socket, NULL, NULL)) == -1) {
                perror("Accepter accept() Error\n");
                exit(1);
            }

            fprintf(accepter_log, "%s: >>> accept connection\n", in_out);

            addElementToList(_accepter_lista_buffer, _M_buffer);

            numero_sender_vivi++;

            pthread_create(&sender_, NULL, &runSender, allocParamSender(communication_socket, _M_buffer));
        }

        fprintf(accepter_log, "%s: ---------- Loop End ----------\n", in_out);

    }

    fprintf(accepter_log, "%s: Socket:\n", in_out);

    /*
     * Chiude le due porte in listening
     */

    close(listening_socket);

    fprintf(accepter_log, "%s: > close listen socket\n", in_out);

    fprintf(accepter_log, "%s: Memoria:\n", in_out);

    free(_parametri_accepter);

    fprintf(accepter_log, "%s: > parametri\n", in_out);

    /*
     * fflush e chiusura del logfile
     */

    fprintf(accepter_log, "%s: Close log\n", in_out);

    fflush(accepter_log);

    fclose(accepter_log);

    return NULL;

}

static parametri_dispatcher allocParamDispatcher(buffer _dispatcher_buffer_in, list _dispatcher_lista_buffer_out) {

    parametri_dispatcher _M_parametri_dispatcher = (parametri_dispatcher) malloc(sizeof (struct parametri_dispatcher_t));

    _M_parametri_dispatcher->dispatcher_buffer_in = _dispatcher_buffer_in;
    _M_parametri_dispatcher->dispatcher_lista_buffer_out = _dispatcher_lista_buffer_out;

    return _M_parametri_dispatcher;

}

/*
 * Quando il dispatcher riceve il POISON_ID allora:
 * - lo inoltra nel sender buffer associato e solo in quello
 * - attende che tutti i receiver siano morti
 * - elimina il receiver buffer associato
 * - fa il fflush del file di log e lo chiude
 * - decrementa numero_dispatcher_vivi
 * - termina
 */

static void *runDispatcher(void *_parametri_dispatcher) {

    parametri_dispatcher _parametri = (parametri_dispatcher) _parametri_dispatcher;

    buffer _dispatcher_buffer_in = _parametri->dispatcher_buffer_in;
    list _dispatcher_lista_buffer_out;

    int _dispatcher_id;
    int _sender_buffer_id;

    list_iterator it_out;
    buffer _sender_buffer;

    msg _msg;

    gsl_matrix *_dispatcher_adjacency_matrix;

    FILE *dispatcher_adjacency_matrix_file;

    FILE *dispatcher_log;

    char dispatcher_log_file_name[FILENAME_MAX];

    _dispatcher_adjacency_matrix = gsl_matrix_alloc(6, 6);
    dispatcher_adjacency_matrix_file = fopen("adjacency_matrix.txt", "r");
    gsl_matrix_fscanf(dispatcher_adjacency_matrix_file, _dispatcher_adjacency_matrix);
    fclose(dispatcher_adjacency_matrix_file);

    /*
    gsl_matrix_fprintf(stdout, adjacency_matrix, "%f");
    printf("valore 0,0 = %f\n",gsl_matrix_get(adjacency_matrix,0,0));
    printf("valore 1,0 = %f\n",gsl_matrix_get(adjacency_matrix,1,0));
    printf("valore 0,1 = %f\n",gsl_matrix_get(adjacency_matrix,0,1));
    printf("valore 1,1 = %f\n",gsl_matrix_get(adjacency_matrix,1,1));
     */

    /* Apre il file di log impostando il suo id a quello del buffer associato */

    _dispatcher_id = getBufferID(_dispatcher_buffer_in);

    sprintf(dispatcher_log_file_name, "dispatcher_%d.log", _dispatcher_id);

    dispatcher_log = fopen(dispatcher_log_file_name, "w");

    if (dispatcher_log == NULL) {
        perror("Errore apertura file: \"dispatcher.log\"\n");
        exit(1);
    };

    setbuf(dispatcher_log, NULL);

    fprintf(dispatcher_log, "D%d: Open log\n", _dispatcher_id);

    /*sleep(1);*/

    fprintf(dispatcher_log, "D%d: Parametri di invocazione:\n", _dispatcher_id);

    _dispatcher_lista_buffer_out = _parametri->dispatcher_lista_buffer_out;

    fprintf(dispatcher_log, "D%d: Wait:\n", _dispatcher_id);

    fprintf(dispatcher_log, "D%d: > getFromBuffer(receiver_buffer)\n", _dispatcher_id);

    /*
     * Riceve un messaggio e lo pone in tutti i sender buffer finchè non riceve il POISON_ID
     */

    _msg = getFromBuffer_blocking(_dispatcher_buffer_in);

    fprintf(dispatcher_log, "D%d: ----------- Loop Start -----------\n", _dispatcher_id);

    while (_msg->msg_id != POISON_ID) {

        it_out = allocListIterator(_dispatcher_lista_buffer_out);

        while (hasNextList(it_out) == 1) {

            int i;

            _sender_buffer = nextElementFromList(it_out);

            _sender_buffer_id = getBufferID(_sender_buffer);

            i = gsl_matrix_get(_dispatcher_adjacency_matrix, _sender_buffer_id - 1, _dispatcher_id - 1);

            /* i = 1;*/

            /* fprintf(dispatcher_log, "sender %d receiver %d matrix_get %d\n", _sender_buffer_id, _dispatcher_id, i); */

            if (i == 1) {

                putInBuffer_blocking(_sender_buffer, _msg);

                fprintf(dispatcher_log, "D%d: > (%d,%d,%f,%f) > %d\n", _dispatcher_id, _msg->msg_id, _msg->msg_time, _msg->msg_x, _msg->msg_z, _sender_buffer_id);
            }
        }

        freeListIterator(it_out);

        fprintf(dispatcher_log, "D%d: ----------------------------------\n", _dispatcher_id);

        _msg = getFromBuffer_blocking(_dispatcher_buffer_in);
    }

    fprintf(dispatcher_log, "D%d: ------------ Loop End ------------\n", _dispatcher_id);

    /*
     * Inoltra nel sender buffer associato e solo in quello il POISON_ID
     */

    it_out = allocListIterator(_dispatcher_lista_buffer_out);

    _sender_buffer = nextElementFromList(it_out);

    _sender_buffer_id = getBufferID(_sender_buffer);

    while (_sender_buffer_id != _dispatcher_id) {

        _sender_buffer = nextElementFromList(it_out);

        _sender_buffer_id = getBufferID(_sender_buffer);
    }

    putInBuffer_blocking(_sender_buffer, _msg);

    fprintf(dispatcher_log, "D%d: Signal:\n", _dispatcher_id);

    fprintf(dispatcher_log, "D%d: > snd (%d,%d,%f,%f) > %d\n", _dispatcher_id, _msg->msg_id, _msg->msg_time, _msg->msg_x, _msg->msg_z, _sender_buffer_id);

    fprintf(dispatcher_log, "D%d: Wait:\n", _dispatcher_id);

    fprintf(dispatcher_log, "D%d: > numero_receiver_vivi > 0\n", _dispatcher_id);

    /*
     * Attende che tutti i receiver siano morti
     */

    while (numero_receiver_vivi > 0);

    fprintf(dispatcher_log, "D%d: Memoria:\n", _dispatcher_id);

    /*
     * free delle variabili usate
     */

    freeBuffer(_dispatcher_buffer_in);

    fprintf(dispatcher_log, "D%d: > free(receiver_buffer)\n", _dispatcher_id);

    free(_parametri_dispatcher);

    fprintf(dispatcher_log, "D%d: > parametri\n", _dispatcher_id);

    /*
     * Fa il fflush del file di log e lo chiude
     */

    fprintf(dispatcher_log, "D%d: Close log\n", _dispatcher_id);

    fflush(dispatcher_log);

    fclose(dispatcher_log);

    /*
     * Decrementa numero_dispatcher_vivi e termina
     */

    numero_dispatcher_vivi--;

    return NULL;

}

parametri_router allocParamRouter(char *_router_ip_to_listen, int _port_client_to_router, int _port_router_to_client) {

    parametri_router _M_parametri_router = (struct parametri_router_t *) malloc(sizeof (struct parametri_router_t));

    _M_parametri_router->router_ip_to_listen = _router_ip_to_listen;
    _M_parametri_router->port_router_to_client = _port_router_to_client;
    _M_parametri_router->port_client_to_router = _port_client_to_router;

    return _M_parametri_router;
}

void *runRouter(void *_parametri_router) {

    parametri_router _parametri = (parametri_router) _parametri_router;

    char *_router_ip_to_listen;
    int _port_client_to_router;
    int _port_router_to_client;

    int _router_agents_number;

    list lista_buffer_in_M_;
    list_iterator it_in;
    buffer _buffer_in;

    list lista_buffer_out_M_;

    pthread_t dispatcher_T_;
    pthread_t accepter_in_T_;
    pthread_t accepter_out_T_;

    FILE *router_log;

    /*
     * Apre il file di log
     */

    router_log = fopen("router.log", "w");

    if (router_log == NULL) {
        perror("Errore apertura file: \"router.log\"\n");
        exit(1);
    };

    setbuf(router_log, NULL);

    fprintf(router_log, ": Open log\n");

    fprintf(router_log, ": Parametri di invocazione:\n");

    _router_ip_to_listen = _parametri->router_ip_to_listen;

    fprintf(router_log, ": > _router_ip_to_listen = %s\n", _router_ip_to_listen);

    _port_client_to_router = _parametri->port_client_to_router;

    fprintf(router_log, ": > _port_client_to_router = %d\n", _port_client_to_router);

    _port_router_to_client = _parametri->port_router_to_client;

    fprintf(router_log, ": > _port_router_to_client = %d\n", _port_router_to_client);

    fprintf(router_log, ": Memoria:\n");

    _router_agents_number = 6;

    fprintf(router_log, ": > _router_agents_number = %d\n", _router_agents_number);

    lista_buffer_in_M_ = allocList();

    fprintf(router_log, ": > alloc(lista_buffer_in)\n");

    lista_buffer_out_M_ = allocList();

    fprintf(router_log, ": > alloc(lista_buffer_out)\n");

    /*
     * Fa partire i due accepter (threads)
     */

    pthread_create(&accepter_in_T_, NULL, &runAccepter, allocParamAccepter(_router_agents_number, _router_ip_to_listen, _port_client_to_router, lista_buffer_in_M_, 0));

    fprintf(router_log, ": Thread -< Accepter():\n");

    fprintf(router_log, ": > _accepter_agents_number = %d\n", _router_agents_number);

    fprintf(router_log, ": > _accepter_address = %s\n", _router_ip_to_listen);

    fprintf(router_log, ": > _accepter_port = %d\n", _port_client_to_router);

    fprintf(router_log, ": > _accepter_in_out = %d\n", 0);

    pthread_create(&accepter_out_T_, NULL, &runAccepter, allocParamAccepter(_router_agents_number, _router_ip_to_listen, _port_router_to_client, lista_buffer_out_M_, 1));

    fprintf(router_log, ": Thread -< Accepter():\n");

    fprintf(router_log, ": > _accepter_agents_number = %d\n", _router_agents_number);

    fprintf(router_log, ": > _accepter_address = %s\n", _router_ip_to_listen);

    fprintf(router_log, ": > _accepter_port = %d\n", _port_router_to_client);

    fprintf(router_log, ": > _accepter_in_out = %d\n", 1);

    /*
     * Attende che i due accepter abbiano terminato il loro compito
     */

    fprintf(router_log, ": Wait:\n");

    fprintf(router_log, ": > pthread_join(accepter_out)\n");

    pthread_join(accepter_out_T_, NULL);

    fprintf(router_log, ": > pthread_join(accepter_in)\n");

    pthread_join(accepter_in_T_, NULL);

    fprintf(router_log, ": Loop:\n");

    fprintf(router_log, ": > while hasNext(lista_buffer_in)\n");

    /*
     * Associa ad ogni receiver buffer un dispatcher
     */

    it_in = allocListIterator(lista_buffer_in_M_);

    fprintf(router_log, ": ----------- Loop Start -----------\n");

    while (hasNextList(it_in) == 1) {

        _buffer_in = nextElementFromList(it_in);

        numero_dispatcher_vivi++;

        pthread_create(&dispatcher_T_, NULL, &runDispatcher, allocParamDispatcher(_buffer_in, lista_buffer_out_M_));

        fprintf(router_log, ": >>> create dispatcher thread\n");
    }

    freeListIterator(it_in);

    fprintf(router_log, ": ------------ Loop End ------------\n");

    /* Attende che tutti i sender thread creati siano morti, infatti loro sono gli ultimi a morire */

    fprintf(router_log, ": Wait:\n");

    fprintf(router_log, ": > numero_sender_vivi > 0\n");

    while (numero_sender_vivi > 0);

    fprintf(router_log, ": Memoria:\n");

    /*
     * free delle variabili usate
     */

    freeList(lista_buffer_in_M_);
    lista_buffer_in_M_ = NULL;

    fprintf(router_log, ": > free(lista_buffer_in)\n");

    freeList(lista_buffer_out_M_);
    lista_buffer_out_M_ = NULL;

    fprintf(router_log, ": > free(lista_buffer_out)\n");

    free(_parametri_router);

    fprintf(router_log, ": > parametri\n");

    /*
     * fflush e chiusura del logfile
     */

    fprintf(router_log, ": Close log\n");

    fflush(router_log);

    fclose(router_log);

    return NULL;
}

max_stars_repo_name: vittorioc/VirtualAgent

max_issues_repo_name: vittorioc/VirtualAgent

max_forks_repo_name: vittorioc/VirtualAgent


text:
/*
 * fvector.h
 *
 *  Created on: Feb 4, 2014
 *      Author: vbonnici
 */

#ifndef FVECTOR_H_
#define FVECTOR_H_

#include <limits>
#include <vector>
#include "data_ts.h"

#include <gsl/gsl_rng.h>

using namespace dolierlib;

namespace dolierlib{

/*
 * shuffle sequences... remember to swap also their lengths
 */
void shuffle(
		std::vector<dna5_t*> &f_sequences,
		std::vector<usize_t> &f_lengths,
		int ntimes)
{
	srand( (unsigned)time(NULL) );
	size_t f,s;
	//dna5_t *tmp_s; usize_t tmp_l;
	for(int i=0; i<ntimes; i++){
		f = static_cast<size_t>(ceil((rand()/(RAND_MAX +1.0))*(static_cast<double>(f_sequences.size() -1))));
		s = static_cast<size_t>(ceil((rand()/(RAND_MAX +1.0))*(static_cast<double>(f_sequences.size() -1))));
//		std::cout<<f<<"\t"<<s<<"\n";
		std::swap(f_sequences[f], f_sequences[s]);
		std::swap(f_lengths[f], f_lengths[s]);
	}
}

/*
 * swap randomly the first leading_possequences with a forward position
 */
void leading_shuffle(
		std::vector<dna5_t*> &f_sequences,
		std::vector<usize_t> &f_lengths,
		size_t leading_pos)
{
	srand( (unsigned)time(NULL) );

	size_t f;
	for(size_t i=0; i<leading_pos; i++){
		f = static_cast<size_t>(ceil((rand()/(RAND_MAX))*(static_cast<double>(f_sequences.size()- leading_pos -1)))) + leading_pos;
		//f = static_cast<size_t>((gsl_rng_uniform (r) * static_cast<double>(f_sequences.size() - leading_pos))) + leading_pos
//		std::cout<<i<<"\t"<<f<<"\n";

//		print_dna5(f_sequences[i], f_lengths[i]);std::cout<<"\t";print_dna5(f_sequences[f], f_lengths[f]);std::cout<<"\n";

		std::swap(f_sequences[i], f_sequences[f]);
		std::swap(f_lengths[i], f_lengths[f]);

//		print_dna5(f_sequences[i], f_lengths[i]);std::cout<<"\t";print_dna5(f_sequences[f], f_lengths[f]);std::cout<<"\n";
	}

}

/*
 * Get kmers spectra.
 * Only non-zero columns are reported.
 * Let X to be the length of the longest input sequences, then theoretically we have sum_{i=1}^{i<=X}(4^i) features (possible kmers)
 * but some of these theoretical kmers do not appear in the input sequences, so they are excluded from the output spectra.
 * Moreover, if sequences are all of the same length, the working X is X-1.
 *
 */
void fvector(
		std::vector<dna5_t*> &f_sequences,
		std::vector<usize_t> &f_lengths,
		double ***vectors,  //output spectra, it's pointer
		size_t *vector_length, //length of spectra, it's a pointer
		double **weights
		)
{

	DNA5MS_t f_ms(true);
	concat(f_ms, f_sequences, f_lengths, false);
	f_ms.areUndefTeminated = true;

	CsFullyContainer f_ds(true);
		build_fully_ds(f_ms, f_ds);
		//build_strict_ds(f_ms, f_ds);
	Cs5DLIndex f_dlindex(f_ms.seq, f_ms.seq_length, f_ds.SA, f_ms.lengths, f_ms.nof_seqs);
	SASearcher sas(f_ms.seq, f_ds.SA, f_ms.seq_length);

	size_t max_length = 0;
	for(size_t i=0; i<f_lengths.size(); i++){
		if(f_lengths[i] > max_length)
			max_length = f_lengths[i];
	}

	size_t nof_features = 0;
	size_t *nof_kfeatures = new size_t[max_length+1];
	for(size_t i=1; i<=max_length; i++){
		NSAIterator it = NSAIterator::begin(f_ms, f_ds, i);
		nof_kfeatures[i] = nof(it);
		nof_features += nof_kfeatures[i];

		std::cout<<i<<"\t"<<nof_kfeatures[i]<<"\t"<<f_sequences.size()<<"\n";

		//if sequences are all of the same length, the working X is X-1.
		if(nof_kfeatures[i] == f_sequences.size()  &&  i<=max_length){
			if(i>1){
				std::cout<<"#\n";
				max_length = i-1;
				nof_features -= nof_kfeatures[i];
			}
			else{
				max_length = i;
			}
			break;
		}
	}

	//double **table = new double*[f_sequences.size()];
	(*vectors) =  new double*[f_sequences.size()];
	for(size_t i=0; i<f_sequences.size(); i++){
		(*vectors)[i] =  new double[nof_features];
		for(size_t j=0; j<nof_features; j++)
			(*vectors)[i][j] = 0;
	}
	(*weights) =  new double[nof_features];

	size_t c_nf = 0;
	size_t ii;
	for(size_t k=1; k<=max_length; k++){
		NSAIterator it = NSAIterator::begin(f_ms, f_ds, k);
		dna5_t *kmer  = new dna5_t[k];
		while(it.next()){
			it.get_kmer(kmer);
//			std::cout<<c_nf<<"\t";print_dna5(kmer,k); std::cout<<"\n";

			for(ii=it.i_start; ii<it.i_end; ii++){
				(*vectors)[f_dlindex.ss_ids[ii]][c_nf]++;
				(*weights)[c_nf] = k;
			}
			c_nf++;
		}
	}

	*vector_length = nof_features;
}

/*
 * normalize vectors by column, from [0,max] to [min,max]/(max-min)
 */
void normalize(double **matrix, size_t size1, size_t size2)
{
	for(size_t i=0; i<size2; i++){
		double min = std::numeric_limits<double>::max();
		double max = std::numeric_limits<double>::min();

		for(size_t j=0; j<size1; j++){
			if(matrix[j][i]< min)
				min = matrix[j][i];
			if(matrix[j][i] > max)
				max = matrix[j][i];
		}

		double f = max - min;
		for(size_t j=0; j<size1; j++){
			matrix[j][i] = (matrix[j][i] - min) / f;
		}
	}
}

/*
 * normalize vectors by column, from [0,max] to [0,max]/(max)
 */
void normalize_by_sum(double **matrix, size_t size1, size_t size2)
{
	for(size_t i=0; i<size2; i++){
		double sum = 0;
		for(size_t j=0; j<size1; j++){
			sum += matrix[j][i];
		}
		for(size_t j=0; j<size1; j++){
			matrix[j][i] /= sum;
		}
	}
}

/*
 * Keep only those columns having at least one element >= min_value.
 * Unselected columns are putted at the end of the matrix.
 */
void keep_only(
		double **matrix, //input matrix
		size_t size1, //number of rows
		size_t size2, //number of columns
		double min_value, //min value for selection
		size_t *new_size2 //result number of columns
		)
{
	size_t to_size = size2;
	size_t to_pos = 0;
	for(size_t j=0; j<size2; j++){
		bool keep = false;
		for(size_t i=0; i<size1; i++){
			if(matrix[i][j] >= min_value){
				keep = true;
				break;
			}
		}

		if(keep){
			if(j != to_pos){
				for(size_t i=0; i<size1; i++){
					matrix[i][to_pos] = matrix[i][j];
				}
			}
			to_pos++;
		}
		else{
			to_size--;
		}
	}
	*new_size2 = to_size;
}

/*
 * Same behaviour of keep_only, but it swap the weights matrix, too.
 */
void keep_only(double **matrix, size_t size1, size_t size2, double min_value, size_t *new_size2, double *weights)
{
	size_t to_size = size2;
	size_t to_pos = 0;
	for(size_t j=0; j<size2; j++){
		bool keep = false;
		for(size_t i=0; i<size1; i++){
			if(matrix[i][j] >= min_value){
				keep = true;
				break;
			}
		}

		if(keep){
			if(j != to_pos){
				for(size_t i=0; i<size1; i++){
					matrix[i][to_pos] = matrix[i][j];
				}
				weights[to_pos] = weights[j];
			}
			to_pos++;
		}
		else{
			to_size--;
		}
	}
	*new_size2 = to_size;
}

/*
 * Theoretical spectra are vectors if length sum_{i=1}^{i<=X}(4^i).
 * Thus, weights[i] correspond to the specific length i of the releated theoretical kmer.
 *
 */
void get_fvectors_weights(double *weights, size_t length, std::vector<dna5_t*> &f_sequences, std::vector<usize_t> &f_lengths){
//	size_t base = 0;
//	size_t ebase = 0;
//	size_t cbase = 0;
//	for(size_t i=0;i<length; i++){
//		if(cbase == ebase){
//			base++;
//			ebase = static_cast<size_t>(pow(4, base));
//			cbase = 0;
//		}
//		weights[i] = base;
//		cbase++;
//	}

	DNA5MS_t f_ms(true);
	concat(f_ms, f_sequences, f_lengths, false);
	f_ms.areUndefTeminated = true;

	CsFullyContainer f_ds(true);
		build_fully_ds(f_ms, f_ds);
		//build_strict_ds(f_ms, f_ds);
	Cs5DLIndex f_dlindex(f_ms.seq, f_ms.seq_length, f_ds.SA, f_ms.lengths, f_ms.nof_seqs);
	SASearcher sas(f_ms.seq, f_ds.SA, f_ms.seq_length);

	size_t max_length = 0;
	for(size_t i=0; i<f_lengths.size(); i++){
		if(f_lengths[i] > max_length)
			max_length = f_lengths[i];
	}

	size_t nof_features = 0;
	size_t *nof_kfeatures = new size_t[max_length+1];
	for(size_t i=1; i<=max_length; i++){
		NSAIterator it = NSAIterator::begin(f_ms, f_ds, i);
		nof_kfeatures[i] = nof(it);
		nof_features += nof_kfeatures[i];

		std::cout<<i<<"\t"<<nof_kfeatures[i]<<"\t"<<f_sequences.size()<<"\n";

		//if sequences are all of the same length, the working X is X-1.
		if(nof_kfeatures[i] == f_sequences.size()  &&  i<=max_length){
			if(i>1){
				std::cout<<"#\n";
				max_length = i-1;
				nof_features -= nof_kfeatures[i];
			}
			else{
				max_length = i;
			}
			break;
		}
	}

	size_t c_nf = 0;
	size_t ii;
	for(size_t k=1; k<=max_length; k++){
		NSAIterator it = NSAIterator::begin(f_ms, f_ds, k);
		dna5_t *kmer  = new dna5_t[k];
		while(it.next()){
			it.get_kmer(kmer);
//			std::cout<<c_nf<<"\t";print_dna5(kmer,k); std::cout<<"\n";
			for(ii=it.i_start; ii<it.i_end; ii++){
				weights[c_nf] = k;
			}
			c_nf++;
		}
	}
}

void print_matrix(double **matrix, size_t size1, size_t size2){
	for(size_t i=0; i<size1; i++){
		for(size_t j=0; j<size2; j++){
			std::cout<<matrix[i][j]<<" ";
		}
		std::cout<<"\n";
	}
}

}

#endif /* FVECTOR_H_ */

ext: h

max_stars_repo_name: vbonnici/DoLiER

max_issues_repo_name: vbonnici/DoLiER

max_forks_repo_name: vbonnici/DoLiER


text:
/*
* -----------------------------------------------------------------
*  Thermochemistry Library --- thrm_lib.h
*  Version: 1.6180
*  Date: Dec 29, 2010
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  This program is free software: you can redistribute it and/or
*  modify it under the terms of the GNU General Public License as
*  published by the Free Software Foundation, either version 3 of
*  the License, or (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the header file for a library with
*  thermochemistry routines.
* -----------------------------------------------------------------
*/

#ifndef __THRM_LIB_H__
#define __THRM_LIB_H__

#include <gsl/gsl_vector.h>
#include <nvector/nvector_serial.h>

/*
* ------------------------------------------------------------
*  Types: struct thrm_struct, thrm_wrk
* ------------------------------------------------------------
* This structure contains fields of the Chemkin workspace.
*
* last update: Mar 5, 2009
* ------------------------------------------------------------
*/

typedef struct thrm_struct
{
    int    *iwk;   /* integer work vector */
    double *rwk;   /* real work vector */
    char   *cwk;   /* character work vector */
    int leniwk;    /* iwk dimension */
    int lenrwk;    /* rwk dimension */
    int lencwk;    /* cwk dimension */
    int n_e;       /* # of chemical elements */
    int n_s;       /* # of elementary species */
    int n_r;       /* # of elementary reactions */
} thrm_wrk;
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

void thrm_title();

void thrm_mech(int n_e,
               int n_s,
               int n_r,
               char *mech);

int thrm_input(int *steps,
               double *t0,
               double *delta_t,
               double *Tmax,
               double *Tmin,
               double *tol);

thrm_wrk* thrm_alloc();

void thrm_free(void **thrm_bl);

int thrm_init(thrm_wrk *thrm);

void thrm_composition(thrm_wrk *thrm,
                      gsl_vector *phi);

void thrm_page_heading0(FILE *file);

double thrm_h2T(thrm_wrk *thrm,
                double *phi,
                double Tmax,
                double Tmin,
                double tol);

void thrm_temp_meanvar(int Np,
                       gsl_vector **ps,
                       thrm_wrk *thrm,
                       double Tmax,
                       double Tmin,
                       double tol,
                       double *mean,
                       double *var);

double thrm_rrsum(gsl_vector *phi);

double thrm_mfsum(gsl_vector *phi);

int thrm_mfsign(gsl_vector *phi);

int thrm_eqs(double t,
             N_Vector phi,
             N_Vector Rphi,
             void *thrm_data);

int thrm_eqs2(double t,
              N_Vector phi,
              N_Vector Sphi,
              void *thrm_data);

#endif /* __THRM_LIB_H__ */

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <Accelerate/Accelerate.h>

#include "infmcmc.h"

#define nj1 32
#define nk1 32

void infmcmc_initChain(INFCHAIN *C, const int nj, const int nk) {
  const int maxk      = (nk >> 1) + 1;
  const int sspectral = sizeof(fftw_complex) * nj * maxk;
  const int sphysical = sizeof(double      ) * nj * nk;
  //const int obsVecMem = sizeof(double      ) * sizeObsVector;
  FILE *fp;
  unsigned long int seed;

  // Set up variables
  C->nj = nj;
  C->nk = nk;
  //C->sizeObsVector = sizeObsVector;
  C->currentIter = 0;
  C->accepted = 0;
  C->_shortTimeAccProbAvg = 0.0;
  C->_bLow = 0.0;
  C->_bHigh = 1.0;

  // Allocate a ton of memory
  C->currentPhysicalState      = (double *)malloc(sphysical);
  C->avgPhysicalState          = (double *)malloc(sphysical);
  C->varPhysicalState          = (double *)malloc(sphysical);
  C->proposedPhysicalState     = (double *)malloc(sphysical);
  C->_M2                       = (double *)malloc(sphysical);
  //C->currentStateObservations  = (double *)malloc(obsVecMem);
  //C->proposedStateObservations = (double *)malloc(obsVecMem);
  //C->data                      = (double *)malloc(obsVecMem);

  C->currentSpectralState  = (fftw_complex *)fftw_malloc(sspectral);
  C->avgSpectralState      = (fftw_complex *)fftw_malloc(sspectral);
  C->priorDraw             = (fftw_complex *)fftw_malloc(sspectral);
  C->proposedSpectralState = (fftw_complex *)fftw_malloc(sspectral);

  memset(C->currentPhysicalState,  0, sphysical);
  memset(C->avgPhysicalState,      0, sphysical);
  memset(C->varPhysicalState,      0, sphysical);
  memset(C->proposedPhysicalState, 0, sphysical);
  memset(C->_M2,                   0, sphysical);
  memset(C->currentSpectralState,  0, sspectral);
  memset(C->avgSpectralState,      0, sspectral);
  memset(C->proposedSpectralState, 0, sspectral);

  C->accProb = 0.0;
  C->avgAccProb = 0.0;
  C->logLHDCurrentState = 0.0;

  /*
   * Set some default values
   */
  C->alphaPrior = 3.0;
  C->rwmhStepSize = 1e-4;
  C->priorVar = 1.0;
  C->priorStd = 1.0;

  C->r = gsl_rng_alloc(gsl_rng_taus2);

  fp = fopen("/dev/urandom", "rb");

  if (fp != NULL) {
    fread(&seed, sizeof(unsigned long int), 1, fp);
    gsl_rng_set(C->r, seed);
    fclose(fp);
    printf("Using random seed\n");
  }
  else {
    gsl_rng_set(C->r, 0);
    printf("Using zero seed\n");
  }

  C->_c2r = fftw_plan_dft_c2r_2d(nj, nk, C->proposedSpectralState, C->proposedPhysicalState, FFTW_MEASURE);
  C->_r2c = fftw_plan_dft_r2c_2d(nj, nk, C->currentPhysicalState, C->currentSpectralState, FFTW_MEASURE);
}

void infmcmc_freeChain(INFCHAIN *C) {
  // Free all allocated memory used by the chain
  free(C->currentPhysicalState);
  free(C->avgPhysicalState);
  free(C->varPhysicalState);
  free(C->proposedPhysicalState);
  free(C->_M2);
  //free(C->currentStateObservations);
  //free(C->proposedStateObservations);

  fftw_free(C->currentSpectralState);
  fftw_free(C->avgSpectralState);
  fftw_free(C->priorDraw);
  fftw_free(C->proposedSpectralState);

  gsl_rng_free(C->r);
}

void infmcmc_resetChain(INFCHAIN *C) {
  infmcmc_freeChain(C);
  infmcmc_initChain(C, C->nj, C->nk);
}

void infmcmc_writeChain(const INFCHAIN *C, FILE *fp) {
  const int s = C->nj * C->nk;

  fwrite(&(C->nj),                 sizeof(int),    1,         fp);
  fwrite(&(C->nk),                 sizeof(int),    1,         fp);
  fwrite(&(C->currentIter),        sizeof(int),    1,         fp);
  fwrite(C->currentPhysicalState,  sizeof(double), s,         fp);
  fwrite(C->avgPhysicalState,      sizeof(double), s,         fp);
  fwrite(C->varPhysicalState,      sizeof(double), s,         fp);
  fwrite(&(C->logLHDCurrentState), sizeof(double), 1,         fp);
  fwrite(&(C->accProb),            sizeof(double), 1,         fp);
  fwrite(&(C->avgAccProb),         sizeof(double), 1,         fp);
}

void infmcmc_writeChainInfo(const INFCHAIN *C, FILE *fp) {
  fwrite(&(C->nj),             sizeof(int), 1, fp);
  fwrite(&(C->nk),             sizeof(int), 1, fp);
}

void infmcmc_writeVFChain(const INFCHAIN *U, const INFCHAIN *V, FILE *fp) {
  const int s = U->nj * U->nk;

  fwrite(U->currentPhysicalState,  sizeof(double), s,         fp);
  fwrite(U->avgPhysicalState,      sizeof(double), s,         fp);
  fwrite(U->varPhysicalState,      sizeof(double), s,         fp);
  fwrite(V->currentPhysicalState,  sizeof(double), s,         fp);
  fwrite(V->avgPhysicalState,      sizeof(double), s,         fp);
  fwrite(V->varPhysicalState,      sizeof(double), s,         fp);
  fwrite(&(U->logLHDCurrentState), sizeof(double), 1,         fp);
  fwrite(&(U->accProb),            sizeof(double), 1,         fp);
  fwrite(&(U->avgAccProb),         sizeof(double), 1,         fp);
}

void infmcmc_printChain(INFCHAIN *C) {
  printf("Iteration %d\n", C->currentIter);
  printf("-- Length is         %d x %d\n", C->nj, C->nk);
  printf("-- llhd val is       %lf\n", C->logLHDCurrentState);
  printf("-- Acc. prob is      %.10lf\n", C->accProb);
  printf("-- Avg. acc. prob is %.10lf\n", C->avgAccProb);
  printf("-- Beta is           %.10lf\n\n", C->rwmhStepSize);
  //finmcmc_printCurrentState(C);
  //finmcmc_printAvgState(C);
  //finmcmc_printVarState(C);
}

void randomPriorDraw(INFCHAIN *C) {
  int j, k;
  const int maxk = (C->nk >> 1) + 1;
  const int nko2 = C->nk >> 1;
  const int njo2 = C->nj >> 1;
  double xrand, yrand, c;
  //const double one = 1.0;

  c = 4.0 * M_PI * M_PI;

  for(j = 0; j < C->nj; j++) {
    for(k = 0; k < maxk; k++) {
      xrand = gsl_ran_gaussian_ziggurat(C->r, C->priorStd);
      if((j == 0) && (k == 0)) {
        C->priorDraw[0] = 0.0;
      }
      else if((j == njo2) && (k == nko2)) {
        C->priorDraw[maxk*njo2+nko2] = /*C->nj */ xrand / pow(c * ((j * j) + (k * k)), (double)C->alphaPrior/2.0);
      }
      else if((j == 0) && (k == nko2)) {
        C->priorDraw[nko2] = /*C->nj */ xrand / pow(c * ((j * j) + (k * k)), (double)C->alphaPrior/2.0);
      }
      else if((j == njo2) && (k == 0)) {
        C->priorDraw[maxk*njo2] = /*C->nj */ xrand / pow(c * ((j * j) + (k * k)), (double)C->alphaPrior/2.0);
      }
      else {
        xrand /= sqrt(2.0);
        yrand = gsl_ran_gaussian_ziggurat(C->r, C->priorStd) / sqrt(2.0);
        C->priorDraw[maxk*j+k] = /*C->nj */ (xrand + I * yrand) / pow(c * ((j * j) + (k * k)), (double)C->alphaPrior/2.0);
        if(j > njo2) {
          C->priorDraw[maxk*j+k] = conj(C->priorDraw[maxk*(C->nj-j)+k]);
        }
      }
    }
  }
}

void randomDivFreePriorDraw(INFCHAIN *C1, INFCHAIN *C2) {
  int j, k;
  const int maxk = (C1->nk >> 1) + 1;
  const int njo2 = C2->nj >> 1;
  const int nko2 = C1->nk >> 1;
  double modk;

  randomPriorDraw(C1);

  for (j = 0; j < C1->nj; j++) {
    for (k = 0; k < maxk; k++) {
      if (j == 0 && k == 0) {
        C1->priorDraw[0] = 0.0;
        C2->priorDraw[0] = 0.0;
        continue;
      }

      modk = sqrt(j * j + k * k);
      if (j < njo2) {
        C2->priorDraw[maxk*j+k] = -j * C1->priorDraw[maxk*j+k] / modk;
      }
      else if (j > njo2) {
        C2->priorDraw[maxk*j+k] = -(j - C1->nj) * C1->priorDraw[maxk*j+k] / modk;
      }
      else {
        C2->priorDraw[maxk*j+k] = 0.0;
      }

      if (k < nko2) {
        C1->priorDraw[maxk*j+k] *= k / modk;
      }
      else {
        C1->priorDraw[maxk*j+k] = 0.0;
      }
    }
  }
}

void infmcmc_seedWithPriorDraw(INFCHAIN *C) {
  const int size = sizeof(fftw_complex) * C->nj * ((C->nk >> 1) + 1);
  fftw_complex *uk = (fftw_complex *)fftw_malloc(size);
  const fftw_plan p = fftw_plan_dft_c2r_2d(C->nj, C->nk, uk, C->currentPhysicalState, FFTW_ESTIMATE);

  randomPriorDraw(C);
  memcpy(uk, C->priorDraw, size);

  //fixme: put into current spectral state

  fftw_execute(p);
  fftw_destroy_plan(p);
  fftw_free(uk);
}

void infmcmc_seedWithDivFreePriorDraw(INFCHAIN *C1, INFCHAIN *C2) {
  const int size = sizeof(fftw_complex) * C1->nj * ((C1->nk >> 1) + 1);
  fftw_complex *uk = (fftw_complex *)fftw_malloc(size);
  const fftw_plan p = fftw_plan_dft_c2r_2d(C1->nj, C1->nk, uk, C1->currentPhysicalState, FFTW_ESTIMATE);

  randomDivFreePriorDraw(C1, C2);

  memcpy(uk, C1->priorDraw, size);
  fftw_execute_dft_c2r(p, uk, C1->currentPhysicalState);

  memcpy(uk, C2->priorDraw, size);
  fftw_execute_dft_c2r(p, uk, C2->currentPhysicalState);

  memcpy(C1->currentSpectralState, C1->priorDraw, size);
  memcpy(C2->currentSpectralState, C2->priorDraw, size);

  fftw_destroy_plan(p);
  fftw_free(uk);
}

void infmcmc_proposeRWMH(INFCHAIN *C) {
  int j, k;
  const int maxk = (C->nk >> 1) + 1;
  const int N = C->nj * C->nk;
  const double sqrtOneMinusBeta2 = sqrt(1.0 - C->rwmhStepSize * C->rwmhStepSize);
  double *u = (double *)malloc(sizeof(double) * C->nj * C->nk);
  fftw_complex *uk = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * C->nj * maxk);

  memcpy(u, C->currentPhysicalState, sizeof(double) * C->nj * C->nk);
  fftw_execute_dft_r2c(C->_r2c, u, C->currentSpectralState);

  // Draw from prior distribution
  randomPriorDraw(C);

  for(j = 0; j < C->nj; j++) {
    for(k = 0; k < maxk; k++) {
      C->proposedSpectralState[maxk*j+k] = sqrtOneMinusBeta2 * C->currentSpectralState[maxk*j+k] + C->rwmhStepSize * C->priorDraw[maxk*j+k];
      C->proposedSpectralState[maxk*j+k] /= N;
    }
  }

  memcpy(uk, C->proposedSpectralState, sizeof(fftw_complex) * C->nj * maxk);
  fftw_execute_dft_c2r(C->_c2r, uk, C->proposedPhysicalState);
  fftw_free(uk);
  free(u);
}

void infmcmc_adaptRWMHStepSize(INFCHAIN *C, double inc) {
  // Adapt to stay in 20-30% range.
  int adaptFreq = 100;
  double rate;

  if (C->currentIter > 0 && C->currentIter % adaptFreq == 0) {
    rate = (double) C->_shortTimeAccProbAvg / adaptFreq;

    if (rate < 0.2) {
      //C->_bHigh = C->rwmhStepSize;
      //C->rwmhStepSize = (C->_bLow + C->_bHigh) / 2.0;
      C->rwmhStepSize -= inc;
    }
    else if (rate > 0.3) {
      //C->_bLow = C->rwmhStepSize;
      //C->rwmhStepSize = (C->_bLow + C->_bHigh) / 2.0;
      C->rwmhStepSize += inc;
    }

    C->_shortTimeAccProbAvg = 0.0;
  }
  else {
    C->_shortTimeAccProbAvg += C->accProb;
  }
}

void infmcmc_proposeDivFreeRWMH(INFCHAIN *C1, INFCHAIN *C2) {
  int j, k;
  const int maxk = (C1->nk >> 1) + 1;
  const int N = C1->nj * C1->nk;
  const double sqrtOneMinusBeta2 = sqrt(1.0 - C1->rwmhStepSize * C1->rwmhStepSize);
  double *u = (double *)malloc(sizeof(double) * C1->nj * C1->nk);
  fftw_complex *uk = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * C1->nj * maxk);

  memcpy(u, C1->currentPhysicalState, sizeof(double) * C1->nj * C1->nk);
  fftw_execute_dft_r2c(C1->_r2c, u, C1->currentSpectralState);

  memcpy(u, C2->currentPhysicalState, sizeof(double) * C2->nj * C2->nk);
  fftw_execute_dft_r2c(C2->_r2c, u, C2->currentSpectralState);

  // Draw from prior distribution
  randomDivFreePriorDraw(C1, C2);

  for(j = 0; j < C1->nj; j++) {
    for(k = 0; k < maxk; k++) {
      C1->currentSpectralState[maxk*j+k] /= N;
      C2->currentSpectralState[maxk*j+k] /= N;
      C1->proposedSpectralState[maxk*j+k] = sqrtOneMinusBeta2 * C1->currentSpectralState[maxk*j+k] + C1->rwmhStepSize * C1->priorDraw[maxk*j+k];
      //C1->proposedSpectralState[maxk*j+k] /= N;
      C2->proposedSpectralState[maxk*j+k] = sqrtOneMinusBeta2 * C2->currentSpectralState[maxk*j+k] + C2->rwmhStepSize * C2->priorDraw[maxk*j+k];
      //C2->proposedSpectralState[maxk*j+k] /= N;
    }
  }

  memcpy(uk, C1->proposedSpectralState, sizeof(fftw_complex) * C1->nj * maxk);
  fftw_execute_dft_c2r(C1->_c2r, uk, C1->proposedPhysicalState);
  memcpy(uk, C2->proposedSpectralState, sizeof(fftw_complex) * C1->nj * maxk);
  fftw_execute_dft_c2r(C2->_c2r, uk, C2->proposedPhysicalState);

  fftw_free(uk);
  free(u);
}

void infmcmc_updateAvgs(INFCHAIN *C) {
  int j, k;
  double deltar;
  fftw_complex deltac;
  C->currentIter++;

  /*
   Update physical vectors
   */
  if (C->currentIter == 1) {
    for (j = 0; j < C->nj; j++) {
      for (k = 0; k < C->nk; k++) {
        deltar = C->currentPhysicalState[C->nk * j + k] - C->avgPhysicalState[C->nk * j + k];
        C->avgPhysicalState[C->nk * j + k] += (deltar / C->currentIter);
        C->_M2[C->nk * j + k] += (deltar * (C->currentPhysicalState[C->nk * j + k] - C->avgPhysicalState[C->nk * j + k]));
        C->varPhysicalState[C->nk * j + k] = -1.0;
      }
    }
  }
  else {
    for (j = 0; j < C->nj; j++) {
      for (k = 0; k < C->nk; k++) {
        deltar = C->currentPhysicalState[C->nk * j + k] - C->avgPhysicalState[C->nk * j + k];
        C->avgPhysicalState[C->nk * j + k] += (deltar / C->currentIter);
        C->_M2[C->nk * j + k] += (deltar * (C->currentPhysicalState[C->nk * j + k] - C->avgPhysicalState[C->nk * j + k]));
        C->varPhysicalState[C->nk * j + k] = C->_M2[C->nk * j + k] / (C->currentIter - 1);
      }
    }
  }

  /*
   Update spectral vectors
  */
  if (C->currentIter == 1) {
    for (j = 0; j < C->nj; j++) {
      for (k = 0; k < C->nk/2 + 1; k++) {
        deltac = C->currentSpectralState[(C->nk/2 + 1) * j + k] - C->avgSpectralState[(C->nk/2 + 1) * j + k];
        C->avgSpectralState[(C->nk/2 + 1) * j + k] += (deltac / C->currentIter);
      }
    }
  }
  else {
    for (j = 0; j < C->nj; j++) {
      for (k = 0; k < C->nk/2 + 1; k++) {
        deltac = C->currentSpectralState[(C->nk/2 + 1) * j + k] - C->avgSpectralState[(C->nk/2 + 1) * j + k];
        C->avgSpectralState[(C->nk/2 + 1) * j + k] += (deltac / C->currentIter);
      }
    }
  }

  /*
   Update scalars
   */
  C->avgAccProb += ((C->accProb - C->avgAccProb) / C->currentIter);
}

void infmcmc_updateRWMH(INFCHAIN *C, double logLHDOfProposal) {
  double alpha;

  alpha = exp(C->logLHDCurrentState - logLHDOfProposal);

  if (alpha > 1.0) {
    alpha = 1.0;
  }
  // fixme: set acc prob here instead
  if (gsl_rng_uniform(C->r) < alpha) {
    memcpy(C->currentSpectralState, C->proposedSpectralState, sizeof(fftw_complex) * C->nj * ((C->nk >> 1) + 1));
    memcpy(C->currentPhysicalState, C->proposedPhysicalState, sizeof(double) * C->nj * C->nk);
    C->accProb = alpha;
    C->logLHDCurrentState = logLHDOfProposal;
  }

  infmcmc_updateAvgs(C);
}

void infmcmc_updateVectorFieldRWMH(INFCHAIN *C1, INFCHAIN *C2, double logLHDOfProposal) {
  double alpha;

  // log likelihoods will be the same for both chains
	alpha = exp(C1->logLHDCurrentState - logLHDOfProposal);

  if (alpha > 1.0) {
    alpha = 1.0;
  }

  C1->accepted = 0;
  C2->accepted = 0;
  C1->accProb = alpha;
  C2->accProb = alpha;

  if (gsl_rng_uniform(C1->r) < alpha) {
    C1->accepted = 1;
    C2->accepted = 1;
    memcpy(C1->currentSpectralState, C1->proposedSpectralState, sizeof(fftw_complex) * C1->nj * ((C1->nk >> 1) + 1));
    memcpy(C1->currentPhysicalState, C1->proposedPhysicalState, sizeof(double) * C1->nj * C1->nk);
    C1->logLHDCurrentState = logLHDOfProposal;
    memcpy(C2->currentSpectralState, C2->proposedSpectralState, sizeof(fftw_complex) * C2->nj * ((C2->nk >> 1) + 1));
    memcpy(C2->currentPhysicalState, C2->proposedPhysicalState, sizeof(double) * C2->nj * C2->nk);
    C2->logLHDCurrentState = logLHDOfProposal;
  }

  infmcmc_updateAvgs(C1);
  infmcmc_updateAvgs(C2);
}

void infmcmc_setRWMHStepSize(INFCHAIN *C, double beta) {
  C->rwmhStepSize = beta;
}

void infmcmc_setPriorAlpha(INFCHAIN *C, double alpha) {
  C->alphaPrior = alpha;
}

void infmcmc_setPriorVar(INFCHAIN *C, double var) {
  C->priorVar = var;
  C->priorStd = sqrt(var);
}

double L2Field(fftw_complex *uk, int nj, int nk) {
  int j, k;
  const int maxk = (nk >> 1) + 1;
  double sum = 0.0;

  for (j = 0; j < nj; j++) {
    for (k = 0; k < maxk; k++) {
      sum += cabs(uk[maxk*j+k]) * cabs(uk[maxk*j+k]);
    }
  }

  return 2.0 * sum;
}

double infmcmc_L2Current(INFCHAIN *C) {
  return L2Field(C->currentSpectralState, C->nj, C->nk);
}

double infmcmc_L2Proposed(INFCHAIN *C) {
  return L2Field(C->proposedSpectralState, C->nj, C->nk);
}

double infmcmc_L2Prior(INFCHAIN *C) {
  return L2Field(C->priorDraw, C->nj, C->nk);
}

void randomPriorDrawOLD(gsl_rng *r, double PRIOR_ALPHA, fftw_complex *randDrawCoeffs) {
  int j, k;
  double xrand, yrand, c;//, scale;
  //const double one = 1.0;

  c = 4.0 * M_PI * M_PI;

  for(j = 0; j < nj1; j++) {
    for(k = 0; k < nk1/2 + 1; k++) {
      if((j == 0) && (k == 0)) {
        randDrawCoeffs[0] = 0.0;
      }
      else if((j == nj1/2) && (k == nk1/2)) {
        xrand = gsl_ran_gaussian_ziggurat(r, 1.0);
        randDrawCoeffs[(nk1/2 + 1) * nj1/2 + nk1/2] = xrand / pow((c * ((j * j) + (k * k))), (double)PRIOR_ALPHA/2.0);
      }
      else if((j == 0) && (k == nk1/2)) {
        xrand = gsl_ran_gaussian_ziggurat(r, 1.0);
        randDrawCoeffs[nk1/2] = xrand / pow((c * ((j * j) + (k * k))), (double)PRIOR_ALPHA/2.0);
      }
      else if((j == nj1/2) && (k == 0)) {
        xrand = gsl_ran_gaussian_ziggurat(r, 1.0);
        randDrawCoeffs[(nk1/2 + 1) * nj1/2] = xrand / pow((c * ((j * j) + (k * k))), (double)PRIOR_ALPHA/2.0);
      }
      else {
        xrand = gsl_ran_gaussian_ziggurat(r, 1.0) / sqrt(2.0);
        yrand = gsl_ran_gaussian_ziggurat(r, 1.0) / sqrt(2.0);
        randDrawCoeffs[(nk1/2 + 1) * j + k] = (xrand + I * yrand) / pow((c * ((j * j) + (k * k))), (double)PRIOR_ALPHA/2.0);
        if(j > nj1/2) {
          randDrawCoeffs[(nk1/2 + 1) * j + k] = conj(randDrawCoeffs[(nk1/2+1)*(nj1-j)+k]);
        }
      }
    }
  }

  /*
  for(j = 1; j < nj1 / 2; j++) {
    for(k = 0; k < nk1 / 2 + 1; k++) {
      xrand = gsl_ran_gaussian_ziggurat(r, 1.0) / M_SQRT2;
      yrand = gsl_ran_gaussian_ziggurat(r, 1.0) / M_SQRT2;
      scale = pow(c * ((j * j) + (k * k)), (double)PRIOR_ALPHA/2.0);
      randDrawCoeffs[(nk1/2 + 1) * j + k] = (xrand + I * yrand) / scale;
      randDrawCoeffs[(nk1/2 + 1) * (nj1-j) + k] = (xrand - I * yrand) / scale;
    }
  }

  for(k = 1; k < nk1 / 2 + 1; k++) {
    xrand = gsl_ran_gaussian_ziggurat(r, 1.0) / M_SQRT2;
    randDrawCoeffs[k] = (xrand + I * yrand) / pow(c * (k * k), (double)PRIOR_ALPHA/2.0);
  }

  for(k = 0; k < nk1/2; k++) {
    xrand = gsl_ran_gaussian_ziggurat(r, 1.0) / M_SQRT2;
    randDrawCoeffs[(nk1/2+1)*(nj1/2)+k] = (xrand + I * yrand) / pow(c * ((nj1 * nj1 / 4) + (k * k)), (double)PRIOR_ALPHA/2.0);
  }

  randDrawCoeffs[0] = 0.0;

  xrand = gsl_ran_gaussian_ziggurat(r, 1.0);
  randDrawCoeffs[(nk1/2+1)*(nj1/2)+(nk1/2)] = xrand / pow(c * ((j * j) + (k * k)), (double)PRIOR_ALPHA/2.0);
  */
}

void setRWMHStepSize(CHAIN *C, double stepSize) {
  C->rwmhStepSize = stepSize;
}

void resetAvgs(CHAIN *C) {
  /*
    Sets avgPhysicalState to currentPhysicalState and
         avgSpectralState to currentSpectralState
  */
  const size_t size_doublenj = sizeof(double) * C->nj;
  memcpy(C->avgPhysicalState, C->currentPhysicalState, size_doublenj * C->nk);
  memcpy(C->avgSpectralState, C->currentSpectralState, size_doublenj * ((C->nk >> 1) + 1));
}

void resetVar(CHAIN *C) {
  // Sets varPhysicalState to 0
  memset(C->varPhysicalState, 0, sizeof(double) * C->nj * C->nk);
}

void proposeIndependence(CHAIN *C) {
  const int maxk = (C->nk >> 1) + 1;

  memcpy(C->proposedSpectralState, C->priorDraw, sizeof(fftw_complex) * C->nj * maxk);
}

double lsqFunctional(const double * const data, const double * const obsVec, const int obsVecSize, const double obsStdDev) {
  int i;
  double temp1, sum1 = 0.0;

  for(i = 0; i < obsVecSize; i++) {
    temp1 = (data[i] - obsVec[i]);
    sum1 += temp1 * temp1;
  }

  return sum1 / (2.0 * obsStdDev * obsStdDev);
}

void acceptReject(CHAIN *C) {
  const double phi1 = lsqFunctional(C->data, C->currentStateObservations, C->sizeObsVector, C->obsStdDev);
  const double phi2 = lsqFunctional(C->data, C->proposedStateObservations, C->sizeObsVector, C->obsStdDev);
  double tempAccProb = exp(phi1 - phi2);

  if(tempAccProb > 1.0) {
    tempAccProb = 1.0;
  }

  if(gsl_rng_uniform(C->r) < tempAccProb) {
    memcpy(C->currentSpectralState, C->proposedSpectralState, sizeof(fftw_complex) * C->nj * ((C->nk >> 1) + 1));
    memcpy(C->currentPhysicalState, C->proposedPhysicalState, sizeof(double) * C->nj * C->nk);
    memcpy(C->currentStateObservations, C->proposedStateObservations, sizeof(double) * C->sizeObsVector);

    C->accProb = tempAccProb;
    C->currentLSQFunctional = phi2;
  }
  else {
    C->currentLSQFunctional = phi1;
  }
}

//
//void updateChain(CHAIN *C) {
//C->currentIter++;
//void updateAvgs(CHAIN *C);
//void updateVar(CHAIN *C);
/*
int main(void) {
  int nj = 32, nk = 32, sizeObsVector = 0;
  unsigned long int randseed = 0;

  FILE *fp;
  CHAIN *C;

  C = (CHAIN *)malloc(sizeof(CHAIN));

  initChain(C, nj, nk, sizeObsVector, randseed);

  randomPriorDraw2(C);

  fp = fopen("2.dat", "w");
  fwrite(C->priorDraw, sizeof(fftw_complex), nj * (nk / 2 + 1), fp);
  fclose(fp);

  return 0;
}
*/

ext: c

max_stars_repo_name: dmcdougall/mcmclib

max_issues_repo_name: dmcdougall/mcmclib

max_forks_repo_name: dmcdougall/mcmclib


text:
static char help[] =
"Solve nonlinear Liouville-Bratu equation in 2D on a structured-grid.  Option prefix lb_.\n"
"Solves\n"
"  - nabla^2 u - lambda e^u = 0\n"
"on the unit square [0,1]x[0,1] subject to zero Dirichlet boundary conditions.\n"
"Critical value occurs about at lambda = 6.808.  Optional exact solution\n"
"(Liouville 1853) in case lambda=1.0.\n\n";

/* compare:
timer ./bratu2D -snes_monitor -snes_converged_reason -ksp_converged_reason -pc_type mg -da_refine 8

timer ./bratu2D -snes_monitor -snes_converged_reason -ksp_converged_reason -pc_type mg -da_refine 8 -snes_fd_color

timer ./bratu2D -snes_monitor -snes_converged_reason -ksp_converged_reason -pc_type mg -da_refine 8 -snes_mf_operator

timer ./bratu2D -snes_monitor -snes_converged_reason -ksp_converged_reason -pc_type mg -snes_grid_sequence 8

timer ./bratu2D -snes_monitor -snes_converged_reason -lb_showcounts -da_refine 8 -snes_type fas -fas_levels_snes_type ngs -fas_coarse_snes_type ngs

timer mpiexec -n 4 ./bratu2D -snes_monitor -snes_converged_reason -lb_showcounts -da_refine 8 -snes_type fas -fas_levels_snes_type ngs -fas_coarse_snes_type ngs

timer ./bratu2D -snes_monitor -snes_converged_reason -lb_showcounts -da_refine 8 -snes_type fas -fas_levels_snes_type ngs -fas_coarse_snes_type ngs -fas_levels_snes_ngs_sweeps 2

timer ./bratu2D -snes_monitor -snes_converged_reason -lb_showcounts -da_refine 8 -snes_type fas -fas_levels_snes_type ngs -fas_coarse_snes_type newtonls -fas_coarse_ksp_type preonly -fas_coarse_pc_type cholesky

timer ./bratu2D -snes_monitor -snes_converged_reason -lb_showcounts -da_refine 8 -snes_type fas -fas_levels_snes_type ngs -fas_coarse_snes_type newtonls -fas_coarse_ksp_type cg -fas_coarse_pc_type icc -snes_fas_monitor -snes_fas_levels 6
*/

/* excellent evidence of convergence in Liouville exact solution case:
$ for LEV in 3 4 5 6 7 8 9; do ./bratu2D -da_refine $LEV -snes_monitor -snes_fd_color -snes_rtol 1.0e-10 -lb_exact -pc_type mg; done
*/

#include <petsc.h>
#include "../../ch6/poissonfunctions.h"

typedef struct {
    PetscReal lambda;
    PetscBool exact;
    int       residualcount, ngscount;
} BratuCtx;

static PetscReal g_zero(PetscReal x, PetscReal y, PetscReal z, void *ctx) {
    return 0.0;
}

static PetscReal g_liouville(PetscReal x, PetscReal y, PetscReal z, void *ctx) {
    PetscReal r2 = (x + 1.0) * (x + 1.0) + (y + 1.0) * (y + 1.0),
              qq = r2 * r2 + 1.0,
              omega = r2 / (qq * qq);
    return 32.0 * omega;
}

extern PetscErrorCode FormUExact(DMDALocalInfo*, Vec, PoissonCtx*);
extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*, PetscReal **,
                                        PetscReal**, PoissonCtx*);
extern PetscErrorCode NonlinearGS(SNES, Vec, Vec, void*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    DM             da, da_after;
    SNES           snes;
    Vec            u, uexact;
    PoissonCtx     user;
    BratuCtx       bctx;
    DMDALocalInfo  info;
    PetscBool      showcounts = PETSC_FALSE;
    PetscLogDouble flops;
    PetscReal      errinf;

    PetscInitialize(&argc,&argv,NULL,help);
    user.Lx = 1.0;
    user.Ly = 1.0;
    user.Lz = 1.0;
    user.cx = 1.0;
    user.cy = 1.0;
    user.cz = 1.0;
    user.g_bdry = &g_zero;
    bctx.lambda = 1.0;
    bctx.exact = PETSC_FALSE;
    bctx.residualcount = 0;
    bctx.ngscount = 0;
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"lb_","Liouville-Bratu equation solver options",""); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-lambda","coefficient of e^u (reaction) term",
                            "bratu2D.c",bctx.lambda,&(bctx.lambda),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-exact","use case of Liouville exact solution",
                            "bratu2D.c",bctx.exact,&(bctx.exact),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-showcounts","at finish, print numbers of calls to call-back functions",
                            "bratu2D.c",showcounts,&showcounts,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);
    if (bctx.exact) {
        if (bctx.lambda != 1.0) {
            SETERRQ(PETSC_COMM_SELF,1,"Liouville exact solution only implemented for lambda = 1.0\n");
        }
        user.g_bdry = &g_liouville;
    }
    user.addctx = &bctx;

    ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,
                        DMDA_STENCIL_BOX,  // contrast with fish2
                        3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da); CHKERRQ(ierr);
    ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);
    ierr = DMSetFromOptions(da); CHKERRQ(ierr);
    ierr = DMSetUp(da); CHKERRQ(ierr);  // this must be called BEFORE SetUniformCoordinates
    ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetDM(snes,da); CHKERRQ(ierr);
    ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
               (DMDASNESFunction)FormFunctionLocal,&user); CHKERRQ(ierr);
    ierr = SNESSetNGS(snes,NonlinearGS,&user); CHKERRQ(ierr);
    // this is the Jacobian of the Poisson equation, thus ONLY APPROXIMATE
    //     ... consider using -snes_fd_color or -snes_mf_operator
    ierr = DMDASNESSetJacobianLocal(da,
               (DMDASNESJacobian)Poisson2DJacobianLocal,&user); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);

    ierr = DMGetGlobalVector(da,&u); CHKERRQ(ierr);
    ierr = VecSet(u,0.0); CHKERRQ(ierr);  // initialize to zero
    ierr = SNESSolve(snes,NULL,u); CHKERRQ(ierr);
    ierr = DMRestoreGlobalVector(da,&u);CHKERRQ(ierr);
    ierr = DMDestroy(&da); CHKERRQ(ierr);

    if (showcounts) {
        ierr = PetscGetFlops(&flops); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_WORLD,"flops = %.3e,  residual calls = %d,  NGS calls = %d\n",
                           flops,bctx.residualcount,bctx.ngscount); CHKERRQ(ierr);
    }

    ierr = SNESGetDM(snes,&da_after); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da_after,&info); CHKERRQ(ierr);
    if (bctx.exact) {
        ierr = SNESGetSolution(snes,&u); CHKERRQ(ierr);  // SNES owns u; we do not destroy it
        ierr = VecDuplicate(u,&uexact); CHKERRQ(ierr);
        ierr = FormUExact(&info,uexact,&user); CHKERRQ(ierr);
        ierr = VecAXPY(u,-1.0,uexact); CHKERRQ(ierr);    // u <- u + (-1.0) uexact
        ierr = VecDestroy(&uexact); CHKERRQ(ierr);  // no longer needed
        ierr = VecNorm(u,NORM_INFINITY,&errinf); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_WORLD,
                "done on %d x %d grid:   error |u-uexact|_inf = %.3e\n",
                info.mx,info.my,errinf); CHKERRQ(ierr);
    } else {
        ierr = PetscPrintf(PETSC_COMM_WORLD,"done on %d x %d grid ...\n",info.mx,info.my); CHKERRQ(ierr);
    }

    ierr = SNESDestroy(&snes); CHKERRQ(ierr);
    return PetscFinalize();
}

PetscErrorCode FormUExact(DMDALocalInfo *info, Vec u, PoissonCtx* user) {
    PetscErrorCode ierr;
    BratuCtx     *bctx = (BratuCtx*)(user->addctx);
    PetscInt     i, j;
    PetscReal    hx, hy, x, y, **au;
    if (user->g_bdry != &g_liouville) {
        SETERRQ(PETSC_COMM_SELF,1,"exact solution only implemented for g_liouville() boundary conditions\n");
    }
    if (bctx->lambda != 1.0) {
        SETERRQ(PETSC_COMM_SELF,2,"Liouville exact solution only implemented for lambda = 1.0\n");
    }
    hx = 1.0 / (PetscReal)(info->mx - 1);
    hy = 1.0 / (PetscReal)(info->my - 1);
    ierr = DMDAVecGetArray(info->da, u, &au);CHKERRQ(ierr);
    for (j=info->ys; j<info->ys+info->ym; j++) {
        y = j * hy;
        for (i=info->xs; i<info->xs+info->xm; i++) {
            x = i * hx;
            au[j][i] = user->g_bdry(x,y,0.0,bctx);
        }
    }
    ierr = DMDAVecRestoreArray(info->da, u, &au);CHKERRQ(ierr);
    return 0;
}

// compute F(u), the residual of the discretized PDE on the given grid
PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscReal **au,
                                 PetscReal **FF, PoissonCtx *user) {
    PetscErrorCode ierr;
    BratuCtx   *bctx = (BratuCtx*)(user->addctx);
    PetscInt   i, j;
    PetscReal  hx, hy, darea, hxhy, hyhx, x, y;

    hx = 1.0 / (PetscReal)(info->mx - 1);
    hy = 1.0 / (PetscReal)(info->my - 1);
    darea = hx * hy;
    hxhy = hx / hy;
    hyhx = hy / hx;
    for (j = info->ys; j < info->ys + info->ym; j++) {
        y = j * hy;
        for (i = info->xs; i < info->xs + info->xm; i++) {
            if (j==0 || i==0 || i==info->mx-1 || j==info->my-1) {
                x = i * hx;
                FF[j][i] = au[j][i] - user->g_bdry(x,y,0.0,bctx);
            } else {
                FF[j][i] =   hyhx * (2.0 * au[j][i] - au[j][i-1] - au[j][i+1])
                           + hxhy * (2.0 * au[j][i] - au[j-1][i] - au[j+1][i])
                           - darea * bctx->lambda * PetscExpScalar(au[j][i]);
            }
        }
    }
    ierr = PetscLogFlops(12.0 * info->xm * info->ym); CHKERRQ(ierr);
    (bctx->residualcount)++;
    return 0;
}

// do nonlinear Gauss-Seidel (processor-block) sweeps on
//     F(u) = b
PetscErrorCode NonlinearGS(SNES snes, Vec u, Vec b, void *ctx) {
    PetscErrorCode ierr;
    PetscInt       i, j, k, maxits, totalits=0, sweeps, l;
    PetscReal      atol, rtol, stol, hx, hy, darea, hxhy, hyhx, x, y,
                   **au, **ab, bij, uu, phi0, phi, dphidu, s;
    DM             da;
    DMDALocalInfo  info;
    PoissonCtx     *user = (PoissonCtx*)(ctx);
    BratuCtx       *bctx = (BratuCtx*)(user->addctx);
    Vec            uloc;

    ierr = SNESNGSGetSweeps(snes,&sweeps);CHKERRQ(ierr);
    ierr = SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&maxits);CHKERRQ(ierr);
    ierr = SNESGetDM(snes,&da);CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);

    hx = 1.0 / (PetscReal)(info.mx - 1);
    hy = 1.0 / (PetscReal)(info.my - 1);
    darea = hx * hy;
    hxhy = hx / hy;
    hyhx = hy / hx;

    ierr = DMGetLocalVector(da,&uloc);CHKERRQ(ierr);
    for (l=0; l<sweeps; l++) {
        ierr = DMGlobalToLocalBegin(da,u,INSERT_VALUES,uloc);CHKERRQ(ierr);
        ierr = DMGlobalToLocalEnd(da,u,INSERT_VALUES,uloc);CHKERRQ(ierr);
        ierr = DMDAVecGetArray(da,uloc,&au);CHKERRQ(ierr);
        if (b) {
            ierr = DMDAVecGetArrayRead(da,b,&ab); CHKERRQ(ierr);
        }
        for (j = info.ys; j < info.ys + info.ym; j++) {
            y = j * hy;
            for (i = info.xs; i < info.xs + info.xm; i++) {
                if (j==0 || i==0 || i==info.mx-1 || j==info.my-1) {
                    x = i * hx;
                    au[j][i] = user->g_bdry(x,y,0.0,bctx);
                } else {
                    if (b)
                        bij = ab[j][i];
                    else
                        bij = 0.0;
                    // do pointwise Newton iterations on scalar function
                    //   phi(u) =   hyhx * (2 u - au[j][i-1] - au[j][i+1])
                    //            + hxhy * (2 u - au[j-1][i] - au[j+1][i])
                    //            - darea * lambda * e^u - bij
                    uu = au[j][i];
                    phi0 = 0.0;
                    for (k = 0; k < maxits; k++) {
                        phi =   hyhx * (2.0 * uu - au[j][i-1] - au[j][i+1])
                              + hxhy * (2.0 * uu - au[j-1][i] - au[j+1][i])
                              - darea * bctx->lambda * PetscExpScalar(uu) - bij;
                        if (k == 0)
                             phi0 = phi;
                        dphidu = 2.0 * (hyhx + hxhy)
                                 - darea * bctx->lambda * PetscExpScalar(uu);
                        s = - phi / dphidu;     // Newton step
                        uu += s;
                        totalits++;
                        if (   atol > PetscAbsReal(phi)
                            || rtol*PetscAbsReal(phi0) > PetscAbsReal(phi)
                            || stol*PetscAbsReal(uu) > PetscAbsReal(s)    ) {
                            break;
                        }
                    }
                    au[j][i] = uu;
                }
            }
        }
        ierr = DMDAVecRestoreArray(da,uloc,&au);CHKERRQ(ierr);
        ierr = DMLocalToGlobalBegin(da,uloc,INSERT_VALUES,u);CHKERRQ(ierr);
        ierr = DMLocalToGlobalEnd(da,uloc,INSERT_VALUES,u);CHKERRQ(ierr);
    }
    ierr = DMRestoreLocalVector(da,&uloc);CHKERRQ(ierr);
    if (b) {
        ierr = DMDAVecRestoreArrayRead(da,b,&ab);CHKERRQ(ierr);
    }
    ierr = PetscLogFlops(21.0 * totalits); CHKERRQ(ierr);
    (bctx->ngscount)++;
    return 0;
}


text:
#include <petsc.h>
#include <petscmath.h>
#include "compressibleFlow.h"
#include "mesh.h"
#include "petscdmplex.h"
#include "petscts.h"

typedef struct {
    PetscInt dim;
    PetscReal k;
    PetscReal gamma;
    PetscReal Rgas;
    PetscReal L;
} Constants;

typedef struct {
    Constants constants;
    FlowData flowData;
} ProblemSetup;

static PetscReal ComputeTExact( PetscReal time, const PetscReal xyz[], Constants *constants, PetscReal rho){

    PetscReal cv = constants->gamma*constants->Rgas/(constants->gamma - 1) - constants->Rgas;

    PetscReal alpha = constants->k/(rho*cv);
    PetscReal Tinitial = 100.0;
    PetscReal T = 300.0;
    for(PetscReal n =1; n < 2000; n ++){
        PetscReal Bn = -Tinitial*2.0*(-1.0 + PetscPowReal(-1.0, n))/(n*PETSC_PI);
        T += Bn*PetscSinReal(n * PETSC_PI*xyz[0]/constants->L)*PetscExpReal(-n*n*PETSC_PI*PETSC_PI*alpha*time/(PetscSqr(constants->L)));
    }

    return T;
}

static PetscErrorCode InitialConditions(PetscInt dim, PetscReal time, const PetscReal xyz[], PetscInt Nf, PetscScalar *node, void *ctx) {
    PetscFunctionBeginUser;

    Constants *constants = (Constants *)ctx;

    PetscReal T = ComputeTExact(time, xyz, constants, 1.0);

    PetscReal u = 0.0;
    PetscReal v = xyz[0];
    PetscReal rho = 1.0;
    PetscReal p= rho*constants->Rgas*T;
    PetscReal e = p/((constants->gamma - 1.0)*rho);
    PetscReal eT = e + 0.5*(u*u + v*v);

    node[RHO] = rho;
    node[RHOE] = rho*eT;
    node[RHOU + 0] = rho*u;
    node[RHOU + 1] = rho*v;

    PetscFunctionReturn(0);
}

static PetscReal computeTemperature(PetscInt dim, const PetscScalar* conservedValues, PetscReal gamma, PetscReal Rgas){
    PetscReal density = conservedValues[RHO];
    PetscReal totalEnergy = conservedValues[RHOE]/density;

    // Get the velocity in this direction
    PetscReal speedSquare = 0.0;
    for (PetscInt d =0; d < dim; d++){
        speedSquare += PetscSqr(conservedValues[RHOU + d]/density);
    }

    // assumed eos
    PetscReal internalEnergy = (totalEnergy) - 0.5 * speedSquare;
    PetscReal p = (gamma - 1.0)*density*internalEnergy;

    PetscReal T = p/(Rgas*density);

    return T;
}

static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx) {
    PetscFunctionBeginUser;
    PetscErrorCode     ierr;

    FlowData flowData = (FlowData)ctx;

    // Get the DM
    DM dm;
    ierr = TSGetDM(ts, &dm);
    CHKERRQ(ierr);

    ierr = FlowViewFromOptions(flowData, "-sol_view");
    CHKERRQ(ierr);

    ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "TS at %f\n", time);
    CHKERRQ(ierr);

    // Compute the error
    void *exactCtxs[1];
    PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
    PetscDS ds;
    ierr = DMGetDS(dm, &ds);
    CHKERRQ(ierr);

    // Get the exact solution
    ierr = PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &exactCtxs[0]);
    CHKERRQ(ierr);

    // Create an vector to hold the exact solution
    Vec exactVec;
    ierr = VecDuplicate(u, &exactVec);
    CHKERRQ(ierr);
    ierr = DMProjectFunction(dm, time, exactFuncs, exactCtxs, INSERT_ALL_VALUES, exactVec);
    CHKERRQ(ierr);

    ierr = PetscObjectSetName((PetscObject)exactVec, "exact");
    CHKERRQ(ierr);
    ierr = VecViewFromOptions(exactVec, NULL, "-sol_view");
    CHKERRQ(ierr);

    // For each component, compute the l2 norms
    ierr = VecAXPY(exactVec, -1.0, u);
    CHKERRQ(ierr);

    PetscReal ferrors[4];
    ierr = VecSetBlockSize(exactVec, 4);
    CHKERRQ(ierr);

    ierr = PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t\n", (int)step, (double)time);
    CHKERRQ(ierr);
    ierr = VecStrideNormAll(exactVec, NORM_2, ferrors);
    CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "\tL_2 Error: [%2.3g, %2.3g, %2.3g, %2.3g]\n", (double)(ferrors[0]), (double)(ferrors[1]), (double)(ferrors[2]), (double)(ferrors[3]));
    CHKERRQ(ierr);

    // And the infinity error
    ierr = VecStrideNormAll(exactVec, NORM_INFINITY, ferrors);
    CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "\tL_Inf Error: [%2.3g, %2.3g, %2.3g, %2.3g]\n", (double)ferrors[0], (double)ferrors[1], (double)ferrors[2], (double)ferrors[3]);
    CHKERRQ(ierr);

    ierr = PetscObjectSetName((PetscObject)exactVec, "error");
    CHKERRQ(ierr);
    ierr = VecViewFromOptions(exactVec, NULL, "-sol_view");
    CHKERRQ(ierr);

    ierr = VecDestroy(&exactVec);
    CHKERRQ(ierr);
    PetscFunctionReturn(0);
}
//
//static PetscErrorCode FillBoundary(PetscInt dim, PetscInt dof, const PetscFVFaceGeom *faceGeom, const PetscFVCellGeom *cellGeom, const PetscFVCellGeom *cellGeomG, const PetscScalar *a_xI, const PetscScalar *a_xGradI, const PetscScalar *a_xG,  PetscScalar *a_xGradG, void *ctx){
//
//    for(PetscInt pd = 0; pd < dof; ++pd) {
//        PetscReal dPhidS = a_xG[pd] - a_xI[pd];
//
//        // over each direction
//        for (PetscInt dir = 0; dir < dim; dir++) {
//            PetscReal dx = (cellGeomG->centroid[dir] - faceGeom->centroid[dir]);
//
//            // If there is a contribution in this direction
//            if (PetscAbs(dx) > 1E-8) {
//                a_xGradG[pd*dim + dir] = dPhidS / (dx);
//            } else {
//                a_xGradG[pd*dim + dir] = 0.0;
//            }
//        }
//    }
//    return 0;
//}
//
///**
// * this function updates the boundaries with the gradient computed from the boundary cell value
// * @param dm
// * @param auxFvm
// * @param gradLocalVec
// * @return
// */
//static PetscErrorCode DMPlexFillGradientInBoundaryFVM(DM dm, PetscFV auxFvm, Vec localXVec, Vec gradLocalVec){
//    PetscFunctionBeginUser;
//    PetscErrorCode ierr;
//
//    // Get the dmGrad
//    DM dmGrad;
//    ierr = VecGetDM(gradLocalVec, &dmGrad);CHKERRQ(ierr);
//
//    // get the problem
//    PetscDS prob;
//    PetscInt nFields;
//    ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
//    ierr = PetscDSGetNumFields(prob, &nFields);CHKERRQ(ierr);
//    if (nFields > 1){
//        SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Cannot fill DMPlexFillGradientInBoundaryFVM with more than a single field." );
//    }
//    PetscInt field;
//    ierr = PetscDSGetFieldIndex(prob, (PetscObject) auxFvm, &field);CHKERRQ(ierr);
//    PetscInt dof;
//    ierr = PetscDSGetFieldSize(prob, field, &dof);CHKERRQ(ierr);
//
//    // Obtaining local cell ownership
//    PetscInt cellStart, cellEnd, cEndInterior;
//    ierr = DMPlexGetHeightStratum(dm, 0, &cellStart, &cellEnd);CHKERRQ(ierr);
//    ierr = DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);CHKERRQ(ierr);
//    cEndInterior = cEndInterior < 0 ? cellEnd: cEndInterior;
//
//    // Get the fvm face and cell geometry
//    Vec cellGeomVec = NULL;/* vector of structs related to cell geometry*/
//    Vec faceGeomVec = NULL;/* vector of structs related to face geometry*/
//    ierr = DMPlexGetGeometryFVM(dm, &faceGeomVec, &cellGeomVec, NULL);CHKERRQ(ierr);
//
//    // get the dm for each geom type
//    DM dmFaceGeom, dmCellGeom;
//    ierr = VecGetDM(faceGeomVec, &dmFaceGeom);CHKERRQ(ierr);
//    ierr = VecGetDM(cellGeomVec, &dmCellGeom);CHKERRQ(ierr);
//
//    // extract the gradLocalVec
//    PetscScalar * gradLocalArray;
//    ierr = VecGetArray(gradLocalVec, &gradLocalArray);CHKERRQ(ierr);
//
//    const PetscScalar* localArray;
//    ierr = VecGetArrayRead(localXVec, &localArray);CHKERRQ(ierr);
//
//    // get the dim
//    PetscInt dim;
//    ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
//
//    // extract the arrays for the face and cell geom, along with their dm
//    const PetscScalar *cellGeomArray;
//    ierr = VecGetArrayRead(cellGeomVec, &cellGeomArray);CHKERRQ(ierr);
//
//    const PetscScalar *faceGeomArray;
//    ierr = VecGetArrayRead(faceGeomVec, &faceGeomArray);CHKERRQ(ierr);
//
//    // March over each face
//    for (PetscInt cell = cellStart; cell < cEndInterior; ++cell) {
//        // Get the face information
//        PetscInt               numFaces;
//        const PetscInt        *faces;
//        ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr);
//        ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr);
//
//        //March over each face
//        for(PetscInt f =0 ; f < numFaces; f++){
//            PetscBool boundary;
//            ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr);
//
//            // if this is on the boundary
//            if(boundary){
//                // get the boundary cell index
//                const PetscInt        *fcells;
//                ierr  = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr);
//                PetscInt side  = (cell != fcells[0]); /* c is on left=0 or right=1 of face */
//                PetscInt ncell = fcells[!side];    /* the neighbor */
//
//                // get the face geom
//                const PetscFVFaceGeom  *faceGeom;
//                ierr  = DMPlexPointLocalRead(dmFaceGeom, faces[f], faceGeomArray, &faceGeom);CHKERRQ(ierr);
//
//                // get the cell centroid information
//                const PetscFVCellGeom       *cellGeom;
//                const PetscFVCellGeom       *cellGeomGhost;
//                ierr  = DMPlexPointLocalRead(dmCellGeom, cell, cellGeomArray, &cellGeom);CHKERRQ(ierr);
//                ierr  = DMPlexPointLocalRead(dmCellGeom, ncell, cellGeomArray, &cellGeomGhost);CHKERRQ(ierr);
//
//                // Read the local point
//                PetscScalar* boundaryGradCellValues;
//                ierr  = DMPlexPointLocalRef(dmGrad, ncell, gradLocalArray, &boundaryGradCellValues);CHKERRQ(ierr);
//
//                const PetscScalar*  cellGradValues;
//                ierr  = DMPlexPointLocalRead(dmGrad, cell, gradLocalArray, &cellGradValues);CHKERRQ(ierr);
//
//                const PetscScalar* boundaryCellValues;
//                ierr  = DMPlexPointLocalRead(dm, ncell, localArray, &boundaryCellValues);CHKERRQ(ierr);
//                const PetscScalar* cellValues;
//                ierr  = DMPlexPointLocalRead(dm, cell, localArray, &cellValues);CHKERRQ(ierr);
//
//                // compute the gradient for the boundary node and pass in
//                //const PetscReal *cellGeom, const PetscReal *cellGeomG, const PetscScalar *a_xI, const PetscScalar *a_xGradI, PetscScalar *a_xG,  PetscScalar *a_xGradG, void *ctx
//                ierr = FillBoundary(dim, dof, faceGeom, cellGeom, cellGeomGhost, cellValues, cellGradValues, boundaryCellValues, boundaryGradCellValues, NULL);CHKERRQ(ierr);
//            }
//        }
//    }
//
//    ierr = VecRestoreArrayRead(localXVec, &localArray);CHKERRQ(ierr);
//    ierr = VecRestoreArray(gradLocalVec, &gradLocalArray);CHKERRQ(ierr);
//    ierr = VecRestoreArrayRead(cellGeomVec, &cellGeomArray);CHKERRQ(ierr);
//    ierr = VecRestoreArrayRead(faceGeomVec, &faceGeomArray);CHKERRQ(ierr);
//
//    PetscFunctionReturn(0);
//}
//
//static PetscErrorCode DiffusionSource(DM dm, PetscReal time, Vec locXVec, Vec globFVec, void *ctx){
//    // Call the flux calculation
//    PetscErrorCode ierr;
//
//    ProblemSetup *setup = (ProblemSetup *)ctx;
//
//    // call the base rhs function eval
//    ierr = DMPlexTSComputeRHSFunctionFVM(dm, time, locXVec, globFVec, setup->flowData);CHKERRQ(ierr);
//
//    // update the aux fields
//    ierr = UpdateAuxFields(NULL, locXVec, setup->flowData);CHKERRQ(ierr);
//
//    Constants constants = setup->constants;
//
//    // get the fvm field assuming that it is the first
//    PetscFV auxFvm;
//    ierr = DMGetField(setup->flowData->auxDm,0, NULL, (PetscObject*)&auxFvm);CHKERRQ(ierr);
//
//    PetscInt dim;
//    ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
//
//    // Get the locXArray
//    const PetscScalar *locXArray;
//    ierr = VecGetArrayRead(locXVec, &locXArray);CHKERRQ(ierr);
//
//    // Get the fvm face and cell geometry
//    Vec cellGeomVec = NULL;/* vector of structs related to cell geometry*/
//    Vec faceGeomVec = NULL;/* vector of structs related to face geometry*/
//
//    // extract the fvm data
//    ierr = DMPlexGetGeometryFVM(setup->flowData->dm, &faceGeomVec, &cellGeomVec, NULL);CHKERRQ(ierr);
//
//    // Get the needed auxDm
//    DM auxFieldGradDM = NULL; /* dm holding the grad information */
//    ierr = DMPlexGetDataFVM(setup->flowData->auxDm, auxFvm, NULL, NULL, &auxFieldGradDM);CHKERRQ(ierr);
//    if(!auxFieldGradDM){
//        SETERRQ(PetscObjectComm((PetscObject)setup->flowData->auxDm), PETSC_ERR_ARG_WRONGSTATE, "The FVM method for aux variables must support computing gradients.");
//    }
//
//    // get the dm for each geom type
//    DM dmFaceGeom, dmCellGeom;
//    ierr = VecGetDM(faceGeomVec, &dmFaceGeom);CHKERRQ(ierr);
//    ierr = VecGetDM(cellGeomVec, &dmCellGeom);CHKERRQ(ierr);
//
//    // extract the arrays for the face and cell geom, along with their dm
//    const PetscScalar *faceGeomArray, *cellGeomArray;
//    ierr = VecGetArrayRead(cellGeomVec, &cellGeomArray);CHKERRQ(ierr);
//    ierr = VecGetArrayRead(faceGeomVec, &faceGeomArray);CHKERRQ(ierr);
//
//    // Obtaining local cell and face ownership
//    PetscInt faceStart, faceEnd;
//    PetscInt cellStart, cellEnd;
//    ierr = DMPlexGetHeightStratum(dm, 1, &faceStart, &faceEnd);CHKERRQ(ierr);
//    ierr = DMPlexGetHeightStratum(dm, 0, &cellStart, &cellEnd);CHKERRQ(ierr);
//
//    // get the ghost label
//    DMLabel ghostLabel;
//    ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr);
//
//    // extract the localFArray from the locFVec
//    PetscScalar *fa;
//    Vec locFVec;
//    ierr = DMGetLocalVector(dm, &locFVec);CHKERRQ(ierr);
//    ierr = VecZeroEntries(locFVec);CHKERRQ(ierr);
//    ierr = VecGetArray(locFVec, &fa);CHKERRQ(ierr);
//
//    // create a global and local grad vector for the auxField
//    Vec gradGlobalVec, gradLocalVec;
//    ierr = DMCreateGlobalVector(auxFieldGradDM, &gradGlobalVec);CHKERRQ(ierr);
//    ierr = VecSet(gradGlobalVec, NAN);CHKERRQ(ierr);
//
//    // compute the global grad values
//    ierr = DMPlexReconstructGradientsFVM(setup->flowData->auxDm, setup->flowData->auxField, gradGlobalVec);CHKERRQ(ierr);
//
//    // Map to a local grad vector
//    ierr = DMCreateLocalVector(auxFieldGradDM, &gradLocalVec);CHKERRQ(ierr);
//    ierr = DMGlobalToLocalBegin(auxFieldGradDM, gradGlobalVec, INSERT_VALUES, gradLocalVec);CHKERRQ(ierr);
//    ierr = DMGlobalToLocalEnd(auxFieldGradDM, gradGlobalVec, INSERT_VALUES, gradLocalVec);CHKERRQ(ierr);
//
//    // fill the boundary conditions
//    ierr = DMPlexFillGradientInBoundaryFVM(setup->flowData->auxDm, auxFvm, setup->flowData->auxField,  gradLocalVec);CHKERRQ(ierr);
//
//    // access the local vector
//    const PetscScalar *localGradArray;
//    ierr = VecGetArrayRead(gradLocalVec,&localGradArray);CHKERRQ(ierr);
//
//    // march over each face
//    for (PetscInt face = faceStart; face < faceEnd; ++face) {
//        PetscFVFaceGeom       *fg;
//        PetscFVCellGeom       *cgL, *cgR;
//        const PetscScalar           *gradL;
//        const PetscScalar           *gradR;
//
//        // make sure that this is a valid face to check
//        PetscInt  ghost, nsupp, nchild;
//        ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr);
//        ierr = DMPlexGetSupportSize(dm, face, &nsupp);CHKERRQ(ierr);
//        ierr = DMPlexGetTreeChildren(dm, face, &nchild, NULL);CHKERRQ(ierr);
//        if (ghost >= 0 || nsupp > 2 || nchild > 0){
//            continue;// skip this face
//        }
//
//        // get the face geometry
//        ierr = DMPlexPointLocalRead(dmFaceGeom, face, faceGeomArray, &fg);CHKERRQ(ierr);
//
//        // Get the left and right cells for this face
//        const PetscInt        *faceCells;
//        ierr = DMPlexGetSupport(dm, face, &faceCells);CHKERRQ(ierr);
//
//        // get the cell geom for the left and right faces
//        ierr = DMPlexPointLocalRead(dmCellGeom, faceCells[0], cellGeomArray, &cgL);CHKERRQ(ierr);
//        ierr = DMPlexPointLocalRead(dmCellGeom, faceCells[1], cellGeomArray, &cgR);CHKERRQ(ierr);
//
//        // extract the cell grad
//        ierr = DMPlexPointLocalRead(auxFieldGradDM, faceCells[0], localGradArray, &gradL);CHKERRQ(ierr);
//        ierr = DMPlexPointLocalRead(auxFieldGradDM, faceCells[1], localGradArray, &gradR);CHKERRQ(ierr);
//
//        PetscInt f = 0;
//
//        // extract the field values
//        PetscScalar *xL, *xR,
//        ierr = DMPlexPointLocalFieldRead(dm, faceCells[0], f, locXArray, &xL);CHKERRQ(ierr);
//        ierr = DMPlexPointLocalFieldRead(dm, faceCells[1], f, locXArray, &xR);CHKERRQ(ierr);
//
//        // Compute the normal grad
//        PetscReal normalGrad = 0.0;
//        PetscInt dof = 0;
//        for (PetscInt d = 0; d < dim; ++d){
//            normalGrad += fg->normal[d]*0.5*(gradL[dof*dim + d] + gradR[dof*dim + d]);
//        }
//
//        // compute the normal heatFlux
//        PetscReal normalHeatFlux = -constants.k *normalGrad;
//
//        // Add to the source terms of f
//        PetscScalar    *fL = NULL, *fR = NULL;
//        ierr = DMLabelGetValue(ghostLabel,faceCells[0],&ghost);CHKERRQ(ierr);
//        if (ghost <= 0) {ierr = DMPlexPointLocalFieldRef(dm, faceCells[0], f, fa, &fL);CHKERRQ(ierr);}
//        ierr = DMLabelGetValue(ghostLabel,faceCells[1],&ghost);CHKERRQ(ierr);
//        if (ghost <= 0) {ierr = DMPlexPointLocalFieldRef(dm, faceCells[1], f, fa, &fR);CHKERRQ(ierr);}
//
//        if(fL){
//            fL[RHOE] -= normalHeatFlux/cgL->volume;
//        }
//        if(fR){
//            fR[RHOE] += normalHeatFlux/cgR->volume;
//        }
//    }
//
//    // Add the new locFVec to the globFVec
//    ierr = VecRestoreArray(locFVec, &fa);CHKERRQ(ierr);
//    ierr = DMLocalToGlobalBegin(dm, locFVec, INSERT_VALUES, globFVec);CHKERRQ(ierr);
//    ierr = DMLocalToGlobalEnd(dm, locFVec, INSERT_VALUES, globFVec);CHKERRQ(ierr);
//    ierr = DMRestoreLocalVector(dm, &locFVec);CHKERRQ(ierr);
//
//    // restore the arrays
//    ierr = VecRestoreArrayRead(gradLocalVec, &localGradArray);CHKERRQ(ierr);
//    ierr = VecRestoreArrayRead(cellGeomVec, &cellGeomArray);CHKERRQ(ierr);
//    ierr = VecRestoreArrayRead(faceGeomVec, &faceGeomArray);CHKERRQ(ierr);
//    ierr = VecRestoreArrayRead(locXVec, &locXArray);CHKERRQ(ierr);
//
//    // destroy grad vectors
//    ierr = VecDestroy(&gradGlobalVec);CHKERRQ(ierr);
//    ierr = VecDestroy(&gradLocalVec);CHKERRQ(ierr);
//
////    VecView(globFVec, PETSC_VIEWER_STDOUT_WORLD);
//
//    PetscFunctionReturn(0);
//}

static PetscErrorCode PhysicsBoundary_Euler(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx) {
    PetscFunctionBeginUser;
    Constants *constants = (Constants *)ctx;

    // compute the centroid location of the real cell
    // Offset the calc assuming the cells are square
    PetscReal x[3];
    for(PetscInt i =0; i < constants->dim; i++){
//        x[i] = c[i] - n[i]*0.5;//TODO:zero boundary
        x[i] = c[i] + n[i]*0.5;//TODO:uniform grad
    }

    // compute the temperature
    PetscReal Tinside =   ComputeTExact(time, x, constants, 1.0);
    PetscReal boundaryValue = 300.0;

    PetscReal T = boundaryValue;//boundaryValue - (Tinside - boundaryValue);//TODO: fix uniform grad

//    PetscReal T = c[0] < constants->L/2.0 ? 300 : 400;

    PetscReal u = 0.0;
    PetscReal v = 0.0;
    PetscReal rho = 1.0;
    PetscReal p= rho*constants->Rgas*T;
    PetscReal e = p/((constants->gamma - 1.0)*rho);
    PetscReal eT = e + 0.5*(u*u + v*v);

    a_xG[RHO] = rho;
    a_xG[RHOE] = rho*eT;
    a_xG[RHOU + 0] = rho*u;
    a_xG[RHOU + 1] = rho*v;

    PetscFunctionReturn(0);
}

static PetscErrorCode PhysicsBoundary_Mirror(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx) {
    PetscFunctionBeginUser;
    Constants *constants = (Constants *)ctx;

    // Offset the calc assuming the cells are square
    for(PetscInt f =0; f < RHOU + constants->dim; f++){
        a_xG[f] = a_xI[f];
    }
    PetscFunctionReturn(0);
}

int main(int argc, char **argv)
{
    // Setup the problem
    Constants constants;

    // sub sonic
    constants.dim = 2;
    constants.L = 0.1;
    constants.gamma = 1.4;
    constants.k = 0.3;
    constants.Rgas = 1.0;

    PetscErrorCode ierr;
    // create the mesh
    // setup the ts
    DM dm;                 /* problem definition */
    TS ts;                 /* timestepper */

    // initialize petsc and mpi
    PetscInitialize(&argc, &argv, NULL, "HELP");

    // Create a ts
    ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);
    ierr = TSSetProblemType(ts, TS_NONLINEAR);CHKERRQ(ierr);
    ierr = TSSetType(ts, TSEULER);CHKERRQ(ierr);
    ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);

    // Create a mesh
    // hard code the problem setup
    PetscReal start[] = {0.0, 0.0};
    PetscReal end[] = {constants.L, constants.L};
    PetscReal nxHeight = 3;
    PetscInt nx[] = {nxHeight, nxHeight};
    ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, constants.dim, PETSC_FALSE, nx, start, end, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr);

    // Setup the flow data
    FlowData flowData;     /* store some of the flow data*/
    ierr = FlowCreate(&flowData);CHKERRQ(ierr);

    // Combine the flow data
    ProblemSetup problemSetup;
    problemSetup.flowData = flowData;
    problemSetup.constants = constants;

    //Setup
    CompressibleFlow_SetupDiscretization(flowData, &dm);

    // Add in the flow parameters
    PetscScalar params[TOTAL_COMPRESSIBLE_FLOW_PARAMETERS];
    params[CFL] = 0.5;
    params[GAMMA] = constants.gamma;
    params[RGAS] = constants.Rgas;
    params[K] = constants.k;

    // set up the finite volume fluxes
    CompressibleFlow_StartProblemSetup(flowData, TOTAL_COMPRESSIBLE_FLOW_PARAMETERS, params);

    // Add in any boundary conditions
    PetscDS prob;
    ierr = DMGetDS(flowData->dm, &prob);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    const PetscInt idsLeft[]= {2, 4};
    ierr = PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall left", "Face Sets", 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler, NULL, 2, idsLeft, &constants);CHKERRQ(ierr);

    const PetscInt idsTop[]= {1, 3};
    ierr = PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "top/bottom", "Face Sets", 0, 0, NULL, (void (*)(void))PhysicsBoundary_Mirror, NULL, 2, idsTop, &constants);CHKERRQ(ierr);

//    ierr = FlowRegisterAuxField(flowData, "Temperature", "T", 1, FV);CHKERRQ(ierr);

    // Complete the problem setup
    ierr = CompressibleFlow_CompleteProblemSetup(flowData, ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Override the flow calc for now
//    ierr = DMTSSetRHSFunctionLocal(flowData->dm, DiffusionSource, &problemSetup);CHKERRQ(ierr);

    // Name the flow field
    ierr = PetscObjectSetName(((PetscObject)flowData->flowField), "Numerical Solution");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Setup the TS
    ierr = TSSetFromOptions(ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSMonitorSet(ts, MonitorError, problemSetup.flowData, NULL);CHKERRQ(ierr);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // set the initial conditions
    PetscErrorCode     (*func[2]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {InitialConditions};
    void* ctxs[1] ={&constants};
    ierr    = DMProjectFunction(flowData->dm,0.0,func,ctxs,INSERT_ALL_VALUES,flowData->flowField);CHKERRQ(ierr);

    // for the mms, add the exact solution
    ierr = PetscDSSetExactSolution(prob, 0, InitialConditions, &constants);CHKERRQ(ierr);

    // Output the mesh
    ierr = DMViewFromOptions(flowData->dm, NULL, "-dm_view");CHKERRQ(ierr);

    // Compute the end time so it goes around once
    // compute dt
    PetscReal cv = constants.gamma*constants.Rgas/(constants.gamma - 1) - constants.Rgas;
    PetscReal dxMin = constants.L/(10 * PetscPowRealInt(2, 2));

    PetscReal alpha = PetscAbs(constants.k/ (cv * 1.0));
    double dt_conduc = 0.00000625;//.3*PetscSqr(dxMin) / alpha;

    PetscPrintf(PETSC_COMM_WORLD, "dt_conduc: %f\n", dt_conduc);
    TSSetTimeStep(ts, dt_conduc);
//    TSSetMaxTime(ts, 0.001);
    TSSetMaxSteps(ts, 600);

    PetscDSView(prob, PETSC_VIEWER_STDOUT_WORLD);

    DMView(flowData->auxDm, PETSC_VIEWER_STDOUT_WORLD);

    ierr = TSSolve(ts,flowData->flowField);CHKERRQ(ierr);

    FlowDestroy(&flowData);
    TSDestroy(&ts);

    return PetscFinalize();

}

max_stars_repo_name: mmcgurn/MattFlowCases

max_issues_repo_name: mmcgurn/MattFlowCases

max_forks_repo_name: mmcgurn/MattFlowCases