/airfoil/airfoil_op.cpp
C++ | 409 lines | 239 code | 90 blank | 80 comment | 37 complexity | ac046ac2c423397fbe0f4be0d2def64f MD5 | raw file
- //
- // auto-generated by op2.m on 30-May-2011 22:03:11
- //
- /*
- Open source copyright declaration based on BSD open source template:
- http://www.opensource.org/licenses/bsd-license.php
- * Copyright (c) 2009-2011, Mike Giles
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- * * Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * * Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- * * The name of Mike Giles may not be used to endorse or promote products
- * derived from this software without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY Mike Giles ''AS IS'' AND ANY
- * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
- * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
- * DISCLAIMED. IN NO EVENT SHALL Mike Giles BE LIABLE FOR ANY
- * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
- * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
- * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
- //
- // Nonlinear airfoil lift calculation
- //
- // Written by Mike Giles, 2010-2011, based on FORTRAN code
- // by Devendra Ghate and Mike Giles, 2005
- //
- //
- // standard headers
- //
- #include <stdlib.h>
- #include <stdio.h>
- #include <string.h>
- #include <math.h>
- #include <assert.h>
- #include <utility>
- #define __NO_STD_VECTOR //use cl::vector
- #include <CL/cl.h>
- // global constants
- //float gam, gm1, cfl, eps, mach, alpha, qinf[4];
- struct global_constants {
- float gam;
- float gm1;
- float cfl;
- float eps;
- float mach;
- float alpha;
- float qinf[4];
- };
- struct global_constants g_const;
- cl_mem g_const_d;
- #include "op_lib.h"
- //#define DIAGNOSTIC 1
- //
- // OP header file
- //
- //
- // op_par_loop declarations
- //
- void op_par_loop_save_soln(char const *, op_set,
- op_arg,
- op_arg );
- void op_par_loop_adt_calc(char const *, op_set,
- op_arg,
- op_arg,
- op_arg,
- op_arg,
- op_arg,
- op_arg );
- void op_par_loop_res_calc(char const *, op_set,
- op_arg,
- op_arg,
- op_arg,
- op_arg,
- op_arg,
- op_arg,
- op_arg,
- op_arg );
- void op_par_loop_bres_calc(char const *, op_set,
- op_arg,
- op_arg,
- op_arg,
- op_arg,
- op_arg,
- op_arg );
- void op_par_loop_update(char const *, op_set,
- op_arg,
- op_arg,
- op_arg,
- op_arg,
- op_arg );
- // kernel routines for parallel loops
- //
- /*
- #include "save_soln.h"
- #include "adt_calc.h"
- #include "res_calc.h"
- #include "bres_calc.h"
- #include "update.h"
- */
- // main program
- void print_array( float *array, int len, const char *file ) {
- FILE *flog;
- flog = fopen( file, "w" );
- for( int i=0; i<len; ++i ) {
- fprintf( flog, "%f\n", array[i] );
- }
- fclose( flog );
- }
- void dump_array( op_dat dat, const char *file ) {
- op_fetch_data( dat );
- print_array( ( float *) dat->data, dat->set->size, file );
- }
- int main(int argc, char **argv){
- int *becell, *ecell, *bound, *bedge, *edge, *cell;
- float *x, *q, *qold, *adt, *res;
- int nnode,ncell,nedge,nbedge,niter;
- float rms;
- // read in grid
- printf("reading in grid \n");
- FILE *fp;
- if ( (fp = fopen("new_grid.dat","r")) == NULL) {
- printf("can't open file new_grid.dat\n"); exit(-1);
- }
- if (fscanf(fp,"%d %d %d %d \n",&nnode, &ncell, &nedge, &nbedge) != 4) {
- printf("error reading from new_grid.dat\n"); exit(-1);
- }
- cell = (int *) malloc(4*ncell*sizeof(int));
- edge = (int *) malloc(2*nedge*sizeof(int));
- ecell = (int *) malloc(2*nedge*sizeof(int));
- bedge = (int *) malloc(2*nbedge*sizeof(int));
- becell = (int *) malloc( nbedge*sizeof(int));
- bound = (int *) malloc( nbedge*sizeof(int));
- x = (float *) malloc(2*nnode*sizeof(float));
- q = (float *) malloc(4*ncell*sizeof(float));
- qold = (float *) malloc(4*ncell*sizeof(float));
- res = (float *) malloc(4*ncell*sizeof(float));
- adt = (float *) malloc( ncell*sizeof(float));
- for (int n=0; n<nnode; n++) {
- if (fscanf(fp,"%f %f \n",&x[2*n], &x[2*n+1]) != 2) {
- printf("error reading from new_grid.dat\n"); exit(-1);
- }
- }
- for (int n=0; n<ncell; n++) {
- if (fscanf(fp,"%d %d %d %d \n",&cell[4*n ], &cell[4*n+1],
- &cell[4*n+2], &cell[4*n+3]) != 4) {
- printf("error reading from new_grid.dat\n"); exit(-1);
- }
- }
- for (int n=0; n<nedge; n++) {
- if (fscanf(fp,"%d %d %d %d \n",&edge[2*n], &edge[2*n+1],
- &ecell[2*n],&ecell[2*n+1]) != 4) {
- printf("error reading from new_grid.dat\n"); exit(-1);
- }
- }
- for (int n=0; n<nbedge; n++) {
- if (fscanf(fp,"%d %d %d %d \n",&bedge[2*n],&bedge[2*n+1],
- &becell[n], &bound[n]) != 4) {
- printf("error reading from new_grid.dat\n"); exit(-1);
- }
- }
- fclose(fp);
- // set constants and initialise flow field and residual
- printf("initialising flow field \n");
- g_const.gam = 1.4f;
- g_const.gm1 = g_const.gam - 1.0f;
- g_const.cfl = 0.9f;
- g_const.eps = 0.05f;
- g_const.mach = 0.4f;
- g_const.alpha = 3.0f*atan(1.0f)/45.0f;
- float p = 1.0f;
- float r = 1.0f;
- float u = sqrt(g_const.gam*p/r)*g_const.mach;
- float e = p/(r*g_const.gm1) + 0.5f*u*u;
- g_const.qinf[0] = r;
- g_const.qinf[1] = r*u;
- g_const.qinf[2] = 0.0f;
- g_const.qinf[3] = r*e;
- for (int n=0; n<ncell; n++) {
- for (int m=0; m<4; m++) {
- q[4*n+m] = g_const.qinf[m];
- res[4*n+m] = 0.0f;
- }
- }
- // OP initialisation
- printf("OP initialisation\n");
- op_init(argc,argv,2);
- g_const_d = op_allocate_constant( &g_const, sizeof( struct global_constants ) );
- // declare sets, pointers, datasets and global constants
- op_set nodes = op_decl_set(nnode, "nodes");
- op_set edges = op_decl_set(nedge, "edges");
- op_set bedges = op_decl_set(nbedge, "bedges");
- op_set cells = op_decl_set(ncell, "cells");
- op_map pedge = op_decl_map(edges, nodes,2,edge, "pedge");
- op_map pecell = op_decl_map(edges, cells,2,ecell, "pecell");
- op_map pbedge = op_decl_map(bedges,nodes,2,bedge, "pbedge");
- op_map pbecell = op_decl_map(bedges,cells,1,becell,"pbecell");
- op_map pcell = op_decl_map(cells, nodes,4,cell, "pcell");
- op_dat p_bound = op_decl_dat(bedges,1,"int" ,bound,"p_bound");
- op_dat p_x = op_decl_dat(nodes ,2,"float",x ,"p_x");
- op_dat p_q = op_decl_dat(cells ,4,"float",q ,"p_q");
- op_dat p_qold = op_decl_dat(cells ,4,"float",qold ,"p_qold");
- op_dat p_adt = op_decl_dat(cells ,1,"float",adt ,"p_adt");
- op_dat p_res = op_decl_dat(cells ,4,"float",res ,"p_res");
- op_decl_const2("gam",1,"float",&g_const.gam );
- op_decl_const2("gm1",1,"float",&g_const.gm1 );
- op_decl_const2("cfl",1,"float",&g_const.cfl );
- op_decl_const2("eps",1,"float",&g_const.eps );
- op_decl_const2("mach",1,"float",&g_const.mach );
- op_decl_const2("alpha",1,"float",&g_const.alpha);
- op_decl_const2("qinf",4,"float",g_const.qinf );
- op_diagnostic_output();
- // main time-marching loop
- niter = 1000;
- for(int iter=1; iter<=niter; iter++) {
- // save old flow solution
- op_par_loop_save_soln("save_soln", cells,
- op_arg_dat(p_q, -1,OP_ID, 4,"float",OP_READ ),
- op_arg_dat(p_qold,-1,OP_ID, 4,"float",OP_WRITE));
- /* if ( iter == 1 ) {
- dump_array( p_qold, "p_qold" );
- }
- */
- #ifdef DIAGNOSTIC
- if (iter==1) {
- dump_array( p_qold, "p_qold" );
- }
- #endif
- //dump_array( p_qold, "p_qold" );
- //op_fetch_data( p_qold );
- //print_array( ( float *) p_qold->data, 4*p_qold->set->size, "p_qold" );
- // print_array( p_q, "p_qold2" );
- // print_array( p_qold, "p_qold" );
- //assert( p_q->data[0] != 0.0f );
- // predictor/corrector update loop
- for(int k=0; k<2; k++) {
- // calculate area/timstep
- op_par_loop_adt_calc("adt_calc",cells,
- op_arg_dat(p_x, 0,pcell, 2,"float",OP_READ ),
- op_arg_dat(p_x, 1,pcell, 2,"float",OP_READ ),
- op_arg_dat(p_x, 2,pcell, 2,"float",OP_READ ),
- op_arg_dat(p_x, 3,pcell, 2,"float",OP_READ ),
- op_arg_dat(p_q, -1,OP_ID, 4,"float",OP_READ ),
- op_arg_dat(p_adt,-1,OP_ID, 1,"float",OP_WRITE));
-
- #ifdef DIAGNOSTIC
- if (iter==1 && k==0) {
- dump_array( p_adt, "p_adt0" );
- }
- if (iter==1 && k==1) {
- dump_array( p_adt, "p_adt1" );
- }
- #endif
- // calculate flux residual
- op_par_loop_res_calc("res_calc",edges,
- op_arg_dat(p_x, 0,pedge, 2,"float",OP_READ),
- op_arg_dat(p_x, 1,pedge, 2,"float",OP_READ),
- op_arg_dat(p_q, 0,pecell,4,"float",OP_READ),
- op_arg_dat(p_q, 1,pecell,4,"float",OP_READ),
- op_arg_dat(p_adt, 0,pecell,1,"float",OP_READ),
- op_arg_dat(p_adt, 1,pecell,1,"float",OP_READ),
- op_arg_dat(p_res, 0,pecell,4,"float",OP_INC ),
- op_arg_dat(p_res, 1,pecell,4,"float",OP_INC ));
- #ifdef DIAGNOSTIC
- if (iter==1 && k==0) {
- dump_array( p_res, "p_res0" );
- }
- if (iter==1 && k==1) {
- dump_array( p_res, "p_res1" );
- }
- #endif
- op_par_loop_bres_calc("bres_calc",bedges,
- op_arg_dat(p_x, 0,pbedge, 2,"float",OP_READ),
- op_arg_dat(p_x, 1,pbedge, 2,"float",OP_READ),
- op_arg_dat(p_q, 0,pbecell,4,"float",OP_READ),
- op_arg_dat(p_adt, 0,pbecell,1,"float",OP_READ),
- op_arg_dat(p_res, 0,pbecell,4,"float",OP_INC ),
- op_arg_dat(p_bound,-1,OP_ID ,1,"int", OP_READ));
- #ifdef DIAGNOSTIC
- if (iter==1 && k==0) {
- dump_array( p_res, "p_res_a0" );
- }
- if (iter==1 && k==0) {
- dump_array( p_res, "p_res_a1" );
- }
- #endif
- // update flow field
- rms = 0.0;
- op_par_loop_update("update",cells,
- op_arg_dat(p_qold,-1,OP_ID, 4,"float",OP_READ ),
- op_arg_dat(p_q, -1,OP_ID, 4,"float",OP_WRITE),
- op_arg_dat(p_res, -1,OP_ID, 4,"float",OP_RW ),
- op_arg_dat(p_adt, -1,OP_ID, 1,"float",OP_READ ),
- op_arg_gbl(&rms,1,"float",OP_INC));
- }
- #ifdef DIAGNOSTIC
- if (iter==1) {
- dump_array( p_q, "p_q1" );
- }
- #endif
- // print iteration history
- rms = sqrt(rms/(float) ncell);
- if (iter%100 == 0)
- printf(" %d %10.5e \n",iter,rms);
- }
- op_timing_output();
- #ifdef DIAGNOSTIC
- dump_array( p_q, "p_q" );
- #endif
- }