testSingle.c 4.35 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
/*******************************************************************************
 * This file is part of SWIFT.
 * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser 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.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 * 
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <pthread.h>
#include <math.h>
#include <float.h>
#include <limits.h>
#include <fenv.h>

/* Conditional headers. */
#ifdef HAVE_LIBZ
    #include <zlib.h>
#endif

/* Local headers. */
#include "swift.h"

/* Ticks per second on this machine. */
#ifndef CPU_TPS
    #define CPU_TPS 2.67e9
#endif

/* Engine policy flags. */
#ifndef ENGINE_POLICY
    #define ENGINE_POLICY engine_policy_none
#endif


/**
 * @brief Main routine that loads a few particles and generates some output.
 *
 */
 
int main ( int argc , char *argv[] ) {

    int k, N = 100;
    struct part p1, p2;
63
    float x, w, dwdx, r2, dx[3] = { 0.0f , 0.0f , 0.0f }, gradw[3];
64
    
65
66
67
    /* Greeting message */
    printf( "This is %s\n", package_description() );

68
69
70
71
72
    /* Init the particles. */
    for ( k = 0 ; k < 3 ; k++ ) {
        p1.a[k] = 0.0f; p1.v[k] = 0.0f; p1.x[k] = 0.0;
        p2.a[k] = 0.0f; p2.v[k] = 0.0f; p2.x[k] = 0.0;
        }
73
74
75
    p1.v[0] = 100.0f;
    p1.id = 0; p2.id = 1;
    p1.density.wcount = 48.0f; p2.density.wcount = 48.0f;
76
77
    p1.rho = 1.0f; p1.mass = 9.7059e-4; p1.h = 0.222871287 / 2;
    p2.rho = 1.0f; p2.mass = 9.7059e-4; p2.h = 0.222871287 / 2;
78
79
    p1.force.c = 0.0040824829f; p1.force.balsara = 0.0f;
    p2.force.c = 58.8972740361f; p2.force.balsara = 0.0f;
80
81
82
83
    p1.u = 1.e-5 / ((const_hydro_gamma - 1.)*p1.rho);
    p2.u = 1.e-5 / ((const_hydro_gamma - 1.)*p2.rho) + 100.0f / ( 33 * p2.mass );
    p1.force.POrho2 = p1.u * ( const_hydro_gamma - 1.0f ) / p1.rho;
    p2.force.POrho2 = p2.u * ( const_hydro_gamma - 1.0f ) / p2.rho;
84
85
    
    /* Dump a header. */
86
87
    printParticle_single( &p1 );
    printParticle_single( &p2 );
88
89
90
91
92
93
    printf( "# r a_1 udt_1 a_2 udt_2\n" );
    
    /* Loop over the different radii. */
    for ( k = 1 ; k <= N ; k++ ) {
    
        /* Set the distance/radius. */
94
        dx[0] = -((float)k)/N * fmaxf( p1.h , p2.h ) * kernel_gamma;
95
96
97
        r2 = dx[0]*dx[0];
        
        /* Clear the particle fields. */
98
99
        p1.a[0] = 0.0f; p1.force.u_dt = 0.0f;
        p2.a[0] = 0.0f; p2.force.u_dt = 0.0f;
100
101
        
        /* Interact the particles. */
102
        runner_iact_force( r2 , dx , p1.h , p2.h , &p1 , &p2 );
103
104
        
        /* Clear the particle fields. */
105
106
        /* p1.rho = 0.0f; p1.density.wcount = 0.0f;
        p2.rho = 0.0f; p2.density.wcount = 0.0f; */
107
108
        
        /* Interact the particles. */
109
        // runner_iact_density( r2 , dx , p1.h , p2.h , &p1 , &p2 );
110
111
112
113
        
        /* Evaluate just the kernel. */
        x = fabsf( dx[0] ) / p1.h;
        kernel_deval( x , &w , &dwdx );
114
115
116
        gradw[0] = dwdx / (p1.h*p1.h*p1.h*p1.h) * dx[0] / sqrtf( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
        gradw[1] = dwdx / (p1.h*p1.h*p1.h*p1.h) * dx[1] / sqrtf( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
        gradw[2] = dwdx / (p1.h*p1.h*p1.h*p1.h) * dx[2] / sqrtf( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
117
118
        
        /* Output the results. */
119
120
121
122
        printf( "%.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e\n" ,
            -dx[0] , p1.a[0] , p1.a[1] , p1.a[2] , p1.force.u_dt , 
            /// -dx[0] , p1.rho , p1.density.wcount , p2.rho , p2.density.wcount ,
            w , dwdx , gradw[0] , gradw[1] , gradw[2] );
123
124
125
126
127
128
129
    
        } /* loop over radii. */
    
    /* All is calm, all is bright. */
    return 0;
    
    }