/home/kapf/tycho_docu/TYCHO.c File Reference

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <time.h>
#include "variables_global.h"
#include "prototypes.h"

Go to the source code of this file.

Functions
int	max_array (int a, int b)
int	main (int argc, char *argv[])

---

Function Documentation

int main	(	int	argc,
		char *	argv[]
	)

This is the main driver of the simulation program. In this function the simulation is initiated, initial conditions are read in or calculated and the main simulation loop is set up. The hydrodynamic simulation is realized with 1D sweeps in all direction and forward/ backward. In addition diffusion and obstacle physics routines are called, when stated in the parameter file.

Definition at line 67 of file TYCHO.c.

References advection, boundary::back, bound, C1_visc, callocate_arrays_global(), cell_volume, counter, counter_restart, dB_Map_ready, diffusion(), dimension, boundary::down, dt, dt_calc(), dt_calc_first(), dt_integrated, endtime_sim, filetype, boundary::front, Gamma, Gamma1, grav_acc, gravity_on_off, hydro_sweeps(), initiate_domain(), initiate_domain_dom(), initiate_domain_marker(), initiate_grid(), intial_soundspeed, boundary::left, make_ic(), make_ics, make_kh_instabilities(), make_sod_ic(), max_array(), max_array_length, minute_in_sec, nmaxx, nmaxy, nmaxz, nminx, nminy, nminz, noise_generator(), number_of_threads, obstacle_heat_conductivity, obstalce_absorption_coefficient, olddt, output_frequency, pre_old_copy(), prepare_the_dB_map(), pressure_on_solid_calc(), pressure_on_solid_reset(), read_dom(), read_ic(), read_marker_file(), read_restart(), read_soundemitter(), read_wind_file(), reset_pressure_integrated(), restart, restart_frequency, boundary::right, S_visc, shock_or_not, small, smallp, smallr, sound_frequency, sound_pressure_level, sound_reflexion, specific_heat_capacity_gas, specific_heat_capacity_obstacle, start_file_reader(), time_sim, boundary::up, velocity_field_analyser(), viscosity_on_off, wind_on_off, with_obstacles, with_one_pulse, with_sound_emitter, write_amira(), write_ic_amira(), write_ic_ifrit(), write_ic_tyc(), write_ic_vtk(), write_ifrit(), write_restart_tyc(), write_tyc(), write_vtk(), x, xmax, xmin, y, ymax, ymin, z, zmax, and zmin.

                                 {

    int tmp1;
    int i, j, k;

#ifndef _OPENMP
    // For performance analysis issues.
    clock_t cputime_sim1, cputime_sim2;
    clock_t cputime_sim1_eta, cputime_sim2_eta;
#endif   

#ifdef _OPENMP
    // For performance analysis issues.
    double cputime_sim1, cputime_sim2;
    double cputime_sim1_eta, cputime_sim2_eta;
#endif   


    specific_heat_capacity_gas = 0.0;
    specific_heat_capacity_obstacle = 0.0;

    double seconds, seconds_eta;

    // File for loading simulation parameter from parameter file.
    char filename[200];

    restart = 0;

    // How many second has a minute ;-)
    minute_in_sec = 60.0;

    //for sound-emitter simulations some initial values
    with_sound_emitter = 0;
    sound_pressure_level = 0.0;
    sound_frequency = 0.0;
    sound_reflexion = 1.0;
    obstalce_absorption_coefficient = 0.0;
    dt_integrated = 0.0;
    dB_Map_ready = 0;
    // Parameterfilename as argument.
    if (argc > 1) {
        sprintf(filename, "%s", argv[1]);
        start_file_reader(filename);
    }

    if (argc == 1) {
        printf("-TYCHO SIMULATION--\n");
        printf("-Linux  Version  --\n");
        printf("-V1.3.0 Oct. 2013--\n");
        printf("-------------------\n");
        printf("Please specify a\n");
        printf("parameter File\n");
        printf("-------------------\n");
        exit(11);
    }

    // restart logic
    if (argc == 3) {
        restart = 1;
        printf("---------------This is a TYCHO Simulation--------------\n");
        printf("---------------------Linux Version --------------------\n");
        printf("---------------Version 1.3.0 released Oct. 2013--------\n");
        printf("The code will restart from last checkpoint\n");
    }

    /*
    Output on the standardout at programme start
     */
    printf("---------------This is a TYCHO Simulation--------------\n");
    printf("---------------------Linux ----------------------------\n");
    printf("---------------Version 1.3.0 released Oct. 2013--------\n");
    printf("The resolution of the simulation is\n%i %i %i\n", x, y, z);
    printf("The domain extend in meter in x,y and z is\n%g %g %g\n", xmax,
            ymax, zmax);
    //to avoid zero time steps
    if ((xmax == 0.0) && (xmax == 0.0) && (zmax == 0.0)) {
        printf("Attention the extend of your computational domain is a singularity!\n");
        printf("We better stop.\n");
        exit(0);
    }
    printf("Your boundary conditions are\n%i %i %i %i %i %i\n",
            bound.down, bound.up, bound.left, bound.right, bound.front, bound.back);
    printf("Your simulation endtime in minutes is\n%g\n", endtime_sim / minute_in_sec);
    if (gravity_on_off == 0) printf("Gravity is off\n");
    if (gravity_on_off == 1) printf("Gravity is on\n");
    if (wind_on_off == 1) printf("A wind marker file will be processed\n");
    if (with_obstacles == 0) printf("No Obstacles within computational domain\n");
    if (with_sound_emitter == 1) printf("With sound emitter\n");
    if (with_one_pulse == 1) printf("With a single sound-shock emitted\n");
    if (advection == 1) printf("A marker field will be advected with the HYDRO solution\n");
    if (viscosity_on_off == 1) printf("Viscosity is switched on\n");
    if (obstacle_heat_conductivity != 0.0) printf("Heat Diffusion is switched on\n");
#ifdef _OPENMP
    printf("You are using %i threads in parallel\n", number_of_threads);
#endif
    printf("-------------------------------------------------------\n");

    //To get output on the standard out
    fflush(stdout);

#ifdef _OPENMP
    // Sets the number of threads specified in the parameter-file.
    omp_set_num_threads(number_of_threads);
#endif

    // Get the maximum number of cells in all three dimension.
    tmp1 = max_array(x, y);
    max_array_length = max_array(tmp1, z);

    // The minimum and maximum index of the cells without ghost cells.
    nminx = 6;
    nmaxx = nminx + x - 1;
    nminy = 6;
    nmaxy = nminy + y - 1;
    nminz = 6;
    nmaxz = nminz + z - 1;

    //Cell Volume
    if (dimension == 3) cell_volume = (xmax / x)*(ymax / y)*(zmax / z);
    //Cell Area
    if (dimension == 2) cell_volume = (xmax / x)*(ymax / y);
    //Cell length
    if (dimension == 1) cell_volume = (xmax / x);

    // Let's start the simulation at the beginning of time.
    time_sim = 0;

    // If Gravitational acceleration the famous constant is defined here.
    if (gravity_on_off == 1) {
        grav_acc = -9.81;
    } else {
        grav_acc = 0.0;
    }

    // Gamma for gas physics.
    Gamma = Gamma1 - 1.0;

    // Small numbers of error handling.
    small = 1.0E-50;
    smallr = 1.0E-50;
    smallp = 1.0E-50;

    //C1_visc and S_visc constant for determine the viscosity with Sutherland's law
    C1_visc = 1.456E-6;
    S_visc = 110.4;

    // The point of origin.
    xmin = 0;
    ymin = 0;
    zmin = 0;

    counter = 0;
    counter_restart = 0;

    //
    intial_soundspeed = 0.0;

    // Dynamic allocation for global variables and arrays.
    callocate_arrays_global(x, y, z, max_array_length);

    //To get output on the standard out
    fflush(stdout);

    // Initiate variables for the hydro solver.
    initiate_grid(x, xmin, xmax, y, ymin, ymax, z, zmin, zmax);

    //To get output on the standard out
    fflush(stdout);

    // Make Initial Conditions or read the from a file in or read only the obstacle domain.
    if (restart == 0) {

        if (advection == 1) {
            read_marker_file(x, y, z);
            initiate_domain_marker(x, y, z);
        }

        if (wind_on_off == 1) {
            read_wind_file(x, y, z);
        }

        if (make_ics == 0) make_ic(x, y, z);
        if (make_ics == 1) read_ic(x, y, z);
        if (make_ics == 2) {
            make_ic(x, y, z);
            read_dom(x, y, z);
        }
        if (make_ics == 3) make_sod_ic(x, y, z);
        if (make_ics == 4) make_kh_instabilities(x, y, z);

        // Initiate variables for the computational domain if obstacles are involved.
        if ((with_obstacles == 1) && (make_ics == 0)) initiate_domain(x, y, z);
        if ((with_obstacles == 1) && (make_ics == 2)) initiate_domain_dom(x, y, z);

        //read in the sound emitter file
        if (with_sound_emitter == 1) read_soundemitter(x, y, z);

        // The initial conditions are written into files.
        if ((filetype == 0) && ((make_ics == 0) || (make_ics == 2) || (make_ics == 3) || (make_ics == 4))) write_ic_tyc(x, y, z);
        if ((filetype == 1) && ((make_ics == 0) || (make_ics == 2) || (make_ics == 3) || (make_ics == 4))) write_ic_vtk(x, y, z);
        if ((filetype == 2) && ((make_ics == 0) || (make_ics == 2) || (make_ics == 3) || (make_ics == 4))) write_ic_amira(x, y, z);
        if ((filetype == 3) && ((make_ics == 0) || (make_ics == 2) || (make_ics == 3) || (make_ics == 4))) write_ic_ifrit(x, y, z);

    }

    if (restart == 1) {

        read_restart(x, y, z);
        printf("Restart files are read.\n");


        //To get output on the standard out
        fflush(stdout);
    }

    // Here the actual hydro solver starts.
    // calculate a first time-step
    dt_calc_first(x, y, z);

    //In the case of a harmonic soundemitter, we have to calculate a new extend, 
    //because of the frequency 100 cells per lambda is quite okay
    //if the extent of the computational domain is large enough this given extend
    //will not be changed
    //if just one soundpulse is emitted (i.e. a weak shock), this is not necessary
    if ((with_sound_emitter == 1) && (with_one_pulse == 0)) {

        double x_max_new, y_max_new, z_max_new;

        x_max_new = xmax;
        y_max_new = ymax;
        z_max_new = zmax;

        double wavelength_divisor = 100.0;
        x_max_new = (intial_soundspeed / sound_frequency) / wavelength_divisor * x;
        if (dimension > 1) y_max_new = (intial_soundspeed / sound_frequency) / wavelength_divisor * y;
        if (dimension > 2) z_max_new = (intial_soundspeed / sound_frequency) / wavelength_divisor * z;

        if ((x_max_new < xmax) || (y_max_new < ymax) || (z_max_new < zmax)) {
            xmax = x_max_new;
            ymax = y_max_new;
            zmax = z_max_new;
            printf("New extend due to harmonic soundsource\nextend in x-direction: %g\nextend in y-direction: %g\nextend in z-direction: %g\n", xmax, ymax, zmax);
            printf("Intial_soundspeed %g [m/s]\n", intial_soundspeed);
            printf("-------------------------------------------------------\n");

            dt_calc_first(x, y, z);
        }
    }

    printf("The beginning time-step: %g [s]\n", dt);
    printf("-------------------------------------------------------\n");

    //To get output on the standard out
    fflush(stdout);

    // For the shock detection we set the variable shock_or_not to zero.
    shock_or_not = 0;

    //the present pressure are copied into the pre_old array
    if ((with_sound_emitter == 1) && (restart == 0)) reset_pressure_integrated(x, y, z);
    if (with_sound_emitter == 1) pre_old_copy(x, y, z);
    if (with_sound_emitter == 1) noise_generator(x, y, z);
    // For performance analysis issues.
    cputime_sim1 = clock();
#ifdef _OPENMP
    cputime_sim1 = omp_get_wtime();
#endif

    // The main Hydro loop, which runs until counter is smaller then a certain integer.
    do {

        // For performance analysis issues.
        cputime_sim1_eta = clock();
#ifdef _OPENMP
        cputime_sim1_eta = omp_get_wtime();
#endif

        //Noise source, i.e. periodic pressure changes
        if (with_sound_emitter == 1) noise_generator(x, y, z);

        //the present pressure are coped into the pre_old array
        if (with_sound_emitter == 1) pre_old_copy(x, y, z);

        // The time-step is now copied to the old time-step.
        olddt = dt;

        // The time-step due to the CFL criteria is determined.
        dt_calc(x, y, z);

        // the arrays where the pressure on the solid is stored 
        // are reseted here
        pressure_on_solid_reset(x, y, z);

        /*
        The Hydro Sweeps forward: In 1D only in x, in 2D in x and y, and
        finally in 3D in the x,y and z direction.
         */
        hydro_sweeps(x, y, z, 0);

        //if heat diffusion is set on
        if (obstacle_heat_conductivity != 0.0) diffusion(x, y, z);

        //Noise source, i.e. periodic pressure changes
        if (with_sound_emitter == 1) noise_generator(x, y, z);

        //the present pressure are coped into the pre_old array
        if (with_sound_emitter == 1) pre_old_copy(x, y, z);
        // The time-step due to the CFL criteria is determined.
        dt_calc(x, y, z);

        // the arrays where the pressure on the solid is stored 
        // are reseted here
        pressure_on_solid_reset(x, y, z);

        /*
         The Hydro Sweeps backward: In 1D only in x, in 2D in x and y, and
         finally in 3D in the x,y and z direction.
         */
        hydro_sweeps(x, y, z, 1);

        // if heat diffusion is set on
        if (obstacle_heat_conductivity != 0.0) diffusion(x, y, z);

        // Velocity Field Analyser
        velocity_field_analyser(x, y, z);

        //Noise source, i.e. periodic pressure changes
        if (with_sound_emitter == 1) noise_generator(x, y, z);

        //the present pressure are coped into the pre_old array
        if (with_sound_emitter == 1) pre_old_copy(x, y, z);

        // Calculate the total pressure on the solid in the computational domain.
        if (with_obstacles == 1) pressure_on_solid_calc(x, y, z);

        // What time is it?
        time_sim = time_sim + 2 * dt;

        // For performance analysis issues.
        cputime_sim2_eta = clock();
#ifdef _OPENMP
        cputime_sim2_eta = omp_get_wtime();
#endif

        seconds_eta = (cputime_sim2_eta - cputime_sim1_eta) / (double) CLOCKS_PER_SEC;
        seconds_eta /= 2.0;

#ifdef _OPENMP
        seconds_eta = (cputime_sim2_eta - cputime_sim1_eta);
        seconds_eta /= 2.0;
#endif
        //check if we have to reset the pressure_integrated array
        //and store the stuff in the dba-map;
        dt_integrated += dt;

        if ((dt_integrated >= 1.0 / sound_frequency) && (with_sound_emitter == 1) && (with_one_pulse == 0)) {

            //make a db map;
            prepare_the_dB_map(x, y, z);
            //reset stuff for the generation of the next dB map
            dt_integrated = 0.0;
            reset_pressure_integrated(x, y, z);
        }

        if ((with_sound_emitter == 1) && (with_one_pulse == 1)) {
            //make a db map;
            prepare_the_dB_map(x, y, z);
        }

        // Write output files at a given frequency
        if (time_sim > output_frequency * counter) {
            if (filetype == 0) write_tyc(x, y, z, counter);
            if (filetype == 1) write_vtk(x, y, z, counter);
            if (filetype == 2) write_amira(x, y, z, counter);
            if (filetype == 3) write_ifrit(x, y, z, counter);

            // For performance analysis issues.
            cputime_sim2 = clock();
#ifdef _OPENMP
            cputime_sim2 = omp_get_wtime();
#endif

            seconds = (cputime_sim2 - cputime_sim1) / (double) CLOCKS_PER_SEC;

#ifdef _OPENMP
            seconds = (cputime_sim2 - cputime_sim1);
#endif

            printf("%3.1f percent done\n", 100.0 / (float) endtime_sim * (float) time_sim);
            printf("time:%g [s], dt:%g [s]\nshock:%i\ncounter:%i\n", time_sim, dt, shock_or_not, counter);

            if (counter > 0) {
                if ((((endtime_sim - time_sim) / dt) * seconds_eta) < 60.0) {
                    printf("Estimated time to simulation end: %4.2f seconds\n", ((endtime_sim - time_sim) / dt) * seconds_eta);
                }
                if (((((endtime_sim - time_sim) / dt) * seconds_eta) > 60.0) && (((((endtime_sim - time_sim) / dt) * seconds_eta) <= 3600.0))) {
                    printf("Estimated time to simulation end: %4.2f minutes\n", ((endtime_sim - time_sim) / dt) * seconds_eta / 60.0);
                }
                if (((((endtime_sim - time_sim) / dt) * seconds_eta) > 3600.0) && ((((endtime_sim - time_sim) / dt) * seconds_eta) <= 3600.0 * 24.0)) {
                    printf("Estimated time to simulation end: %5.2f hours\n", ((endtime_sim - time_sim) / dt) * seconds_eta / 3600.0);
                }
                if ((((endtime_sim - time_sim) / dt) * seconds_eta) > 3600.0 * 24.0) {
                    printf("Estimated time to simulation end: %5.2f days\n", ((endtime_sim - time_sim) / dt) * seconds_eta / (3600.0 * 24.0));
                }
                if (seconds < 60.0) printf("Time since last output: %4.2g [s]\n", seconds);
                if ((seconds > 60.0) && (seconds < 3600)) printf("Time since last output: %4.2g [min]\n", seconds / 60.0);
                if (seconds > 3600.0) printf("Time since last output: %4.2g [h]\n", seconds / 3600.0);
            }
            printf("-------------------------------------------------------\n");

            if (time_sim > restart_frequency * counter_restart) {
                if (filetype != 0) {
                    write_restart_tyc(x, y, z);
                    printf("Restart files are written.\n");
                    counter_restart++;
                }
            }

            counter++;

            //To get output on the standard out
            fflush(stdout);
            // For performance analysis issues.

            cputime_sim1 = clock();
#ifdef _OPENMP
            cputime_sim1 = omp_get_wtime();
#endif
        }


    } while (time_sim <= endtime_sim);

    printf("\n--------------------------");
    printf("\n--------------------------\n");
    printf("Program aborted normally\n");
    printf("End-time of simulation reached\n");
    printf("\n--------------------------");
    printf("\n--------------------------\n");

    return 0;
}

int max_array ( int  a, int  b  )

inline

Definition at line 45 of file TYCHO.c.

Referenced by main().

                                   {
    double tmp;

    if (a < b) tmp = b;
    if (a == b)tmp = b;
    if (a > b) tmp = a;

    return tmp;

}