D
Deleted member 17345
Cryptocurrency Expert
- Joined
- Jan 21, 2022
- Posts
- 3,911
- Reputation
- 2,696
#include <stdio.h>
#include <math.h>
int main(){
/* Grid properties */
const int NX=1000; // Number of x points reduce this to 100 now but change this to 1000 when submiting work
const int NY=1000; // Number of y points
const float xmin=0.0; // Minimum x value
const float xmax=30.0; // Maximum x value
const float ymin=0.0; // Minimum y value
const float ymax=30.0; // Maximum y value
/* Parameters for the Gaussian initial conditions */
const float x0=3.0; // Centre(x)
const float y0=15.0; // Centre
const float sigmax=1.00; // Width(x)
const float sigmay=5.00; // Width
const float sigmax2 = sigmax * sigmax; // Width(x) squared
const float sigmay2 = sigmay * sigmay; // Width squared
/* Boundary conditions */
const float bval_left=0.0; // Left boudnary value
const float bval_right=0.0; // Right boundary value
const float bval_lower=0.0; // Lower boundary
const float bval_upper=0.0; // Upper bounary
/* Time stepping parameters */
const float CFL=0.9; // CFL number
const int nsteps=800; // Number of time steps
/* Velocity */
float velx[NX+2]; // Velocity in x direction will now vary with height so its not a const anymore
const float vely =0.00; // Velocity in y direction
/* Arrays to store variables. These have NX+2 elements
to allow boundary values to be stored at both ends */
float x[NX+2]; // x-axis values
float y[NX+2]; // y-axis values (height z)
float u[NX+2][NY+2]; // Array of u values
float dudt[NX+2][NY+2]; // Rate of change of u
float x2; // x squared (used to calculate iniital conditions)
float y2; // y squared (used to calculate iniital conditions)
/* Calculate distance between points */
float dx = (xmax-xmin) / ( (float) NX);
float dy = (ymax-ymin) / ( (float) NY);
float zZero = 1.0; //paramaters for the shear logarithmic equation
float uStar= 0.2;
float kValue = 0.41;
float velocity_constant = uStar/kValue; //showing this explicity
/* Calculate time step using the CFL condition */
/* The fabs function gives the absolute value in case the velocity is -ve */
float dt[NX+2];
//the two loops can be parallelised without a problem
#pragma omp parallel deafult(none) private(y,velx,dt) //check how to make array private ion open mp
{
/*** Place x points in the middle of the cell ***/
/* LOOP 1 */
for (int i=0; i<NX+2; i++){
x = ( (float) i - 0.5) * dx;
}
/*** Place y points in the middle of the cell ***/
/* LOOP 2 */
for (int j=0; j<NY+2; j++){
y[j] = ( (float) j - 0.5) * dy;
if y[j] > zZero { //will only apply the logarithmic profile once height z is above roughtness length.
velx[j] = log(y[j]/zZero)*velocity_constant;
}else{
velx[j] = 0;
}
dt[j]= CFL / ( (fabs(velx[j]) / dx) + (fabs(vely) / dy)); //this thing will vary with velx changing?
}
}
/*** Report information about the calculation ***/
printf("Grid spacing dx = %g\n", dx);
printf("Grid spacing dy = %g\n", dy);
printf("CFL number = %g\n", CFL);
printf("Time step = %g\n", dt[0]);//chnages with horizontal velocity velocity
printf("No. of time steps = %d\n", nsteps);
printf("End time = %g\n", dt[0]*(float) nsteps);
printf("Distance advected x = %g\n", velx[0]*dt[0]*(float) nsteps); //this will vary as velocity will vary with height
printf("Distance advected y = %g\n", vely*dt*(float) nsteps);
/*** Set up Gaussian initial conditions ***/
/* LOOP 3 */
//here we can use collapse keyword as we have a inner loop thats not related to the outerloop by iteration counts
#pragma omp parallel deafult(none) collapse(2) private(x2,y2) //we have to make x and y2 private to not cause unexpected results due to race conditions
for (int i=0; i<NX+2; i++){
for (int j=0; j<NY+2; j++){
x2 = (x-x0) * (x-x0);
y2 = (y[j]-y0) * (y[j]-y0);
u[j] = exp(-1.0 * ( (x2/(2.0*sigmax2)) + (y2/(2.0*sigmay2))));
}
}
/*** Write array of initial u values out to file ***/
FILE *initialfile;
initialfile = fopen("initial.dat", "w");
/* LOOP 4 */
for (int i=0; i<NX+2; i++){ //cannot paralise because this is data output so it has to be delivered in a specific order an multithreading doesnt gurantee that
for (int j=0; j<NY+2; j++){
fprintf(initialfile, "%g %g %g\n", x, y[j], u[j]);
}
}
fclose(initialfile);
/*** Update solution by looping over time steps ***/
/* LOOP 5 */
for (int m=0; m<nsteps; m++){//cant just use collapse here as some loops inside this loop cannot be paralised
#pragma omp parallel deafult(none)
{
/*** Apply boundary conditions at u[0][:] and u[NX+1][:] ***/
/* LOOP 6 */
for (int j=0; j<NY+2; j++){
u[0][j] = bval_left;
u[NX+1][j] = bval_right;
}
/*** Apply boundary conditions at u[:][0] and u[:][NY+1] ***/
/* LOOP 7 */
for (int i=0; i<NX+2; i++){
u[0] = bval_lower;
u[NY+1] = bval_upper;
}
}
/*** Calculate rate of change of u using leftward difference ***/
/* Loop over points in the domain but not boundary values */
/* LOOP 8 */
//cant parallelise because we have a dependency where we need i-1 elements that could of been executed by other threads causing unexpected results
for (int i=1; i<NX+1; i++){
for (int j=1; j<NY+1; j++){
dudt[j] = -velx[j] * (u[j] - u[i-1][j]) / dx //velx is now variable depending on the height
- vely * (u[j] - u[j-1]) / dy;
}
}
/*** Update u from t to t+dt ***/
/* Loop over points in the domain but not boundary values */
/* LOOP 9 */
#pragma omp parallel for deafult(none) collapse(2) reduction(+ : u)//we have to perform a reduction as it would give innacuare results because here we are accumulating a sum so it has to be done in a specific manner for the result to be valid
for (int i=1; i<NX+1; i++){
for (int j=1; j<NY+1; j++){
u[j] = u[j] + dudt[j] * dt[j];
}
}
} // time loop
/*** Write array of final u values out to file ***/
FILE *finalfile;
finalfile = fopen("final.dat", "w");
/* LOOP 10 */
//cant paralise once again because we need ot deliver data in only one order which isnt guranteed by multithreading
for (int i=0; i<NX+2; i++){
for (int j=0; j<NY+2; j++){
fprintf(finalfile, "%g %g %g\n", x, y[j], u[j]);
}
}
fclose(finalfile);
//creating vertical average we want to plot a average value of u for each x point.
FILE *verticalaverage;
verticalaverage = fopen("VerticalAverage.dat", "w");
/* LOOP 11 */
//cant paralise once again because we need ot deliver data in only one order which isnt guranteed by multithreading
float vertical_avg[NX+2][NX+2][NX+2]; //a array that holds average values of u for every x value
int itemCounter = 0 ;
bool found = false;
// {X VALUE , total U VALUE,total number of u values }
//CONSIDR BOUNDARY THO
for (int i=0; i<NX+2; i++){
for (int j=0; j<NY+2; j++){
//IF FOUND ADD U AND AVG
found = false;
for (int k=0; k<itemCounter; k++){ //check if paricular x value was inserted already
if (vertical_avg[k][0][0]==x){
//add to the array a new value of x
vertical_avg[k][1][0]=vertical_avg[k][1][0] + u[j];
vertical_avg[k][1][1]=vertical_avg[k][1][1] + 1; //increments the count as we added a new value
found =true;
break;//exit loop as its done for this iteration
}
}
if (found == false){
//add new x value into the array
vertical_avg[itemCounter][0][0] = x;
vertical_avg[itemCounter][1][0] = u[j];
vertical_avg[itemCounter][1][1] +=1 ;
itemCounter+=1;
}
}
}
//now we need to divide the total sum by number of elements that created the total value of u to find the average u for each x value.
for (int k= 0 ;k<itermcounter;k++){
vertical_avg[itemCounter][1][0] = vertical_avg[itemCounter][1][0]/ vertical_avg[itemCounter][1][1];
//each field here i now the average value of u for every x value.
fprintf(fverticalaverage, "%g %g %g\n", vertical_avg[itemCounter][0][0], vertical_avg[itemCounter][1][0]);
prints x value on the x axis and the u value of the y axis.
}
fclose(verticalaverage);
return 0;
}
#include <math.h>
int main(){
/* Grid properties */
const int NX=1000; // Number of x points reduce this to 100 now but change this to 1000 when submiting work
const int NY=1000; // Number of y points
const float xmin=0.0; // Minimum x value
const float xmax=30.0; // Maximum x value
const float ymin=0.0; // Minimum y value
const float ymax=30.0; // Maximum y value
/* Parameters for the Gaussian initial conditions */
const float x0=3.0; // Centre(x)
const float y0=15.0; // Centre
const float sigmax=1.00; // Width(x)
const float sigmay=5.00; // Width
const float sigmax2 = sigmax * sigmax; // Width(x) squared
const float sigmay2 = sigmay * sigmay; // Width
squared/* Boundary conditions */
const float bval_left=0.0; // Left boudnary value
const float bval_right=0.0; // Right boundary value
const float bval_lower=0.0; // Lower boundary
const float bval_upper=0.0; // Upper bounary
/* Time stepping parameters */
const float CFL=0.9; // CFL number
const int nsteps=800; // Number of time steps
/* Velocity */
float velx[NX+2]; // Velocity in x direction will now vary with height so its not a const anymore
const float vely =0.00; // Velocity in y direction
/* Arrays to store variables. These have NX+2 elements
to allow boundary values to be stored at both ends */
float x[NX+2]; // x-axis values
float y[NX+2]; // y-axis values (height z)
float u[NX+2][NY+2]; // Array of u values
float dudt[NX+2][NY+2]; // Rate of change of u
float x2; // x squared (used to calculate iniital conditions)
float y2; // y squared (used to calculate iniital conditions)
/* Calculate distance between points */
float dx = (xmax-xmin) / ( (float) NX);
float dy = (ymax-ymin) / ( (float) NY);
float zZero = 1.0; //paramaters for the shear logarithmic equation
float uStar= 0.2;
float kValue = 0.41;
float velocity_constant = uStar/kValue; //showing this explicity
/* Calculate time step using the CFL condition */
/* The fabs function gives the absolute value in case the velocity is -ve */
float dt[NX+2];
//the two loops can be parallelised without a problem
#pragma omp parallel deafult(none) private(y,velx,dt) //check how to make array private ion open mp
{
/*** Place x points in the middle of the cell ***/
/* LOOP 1 */
for (int i=0; i<NX+2; i++){
x = ( (float) i - 0.5) * dx;
}
/*** Place y points in the middle of the cell ***/
/* LOOP 2 */
for (int j=0; j<NY+2; j++){
y[j] = ( (float) j - 0.5) * dy;
if y[j] > zZero { //will only apply the logarithmic profile once height z is above roughtness length.
velx[j] = log(y[j]/zZero)*velocity_constant;
}else{
velx[j] = 0;
}
dt[j]= CFL / ( (fabs(velx[j]) / dx) + (fabs(vely) / dy)); //this thing will vary with velx changing?
}
}
/*** Report information about the calculation ***/
printf("Grid spacing dx = %g\n", dx);
printf("Grid spacing dy = %g\n", dy);
printf("CFL number = %g\n", CFL);
printf("Time step = %g\n", dt[0]);//chnages with horizontal velocity velocity
printf("No. of time steps = %d\n", nsteps);
printf("End time = %g\n", dt[0]*(float) nsteps);
printf("Distance advected x = %g\n", velx[0]*dt[0]*(float) nsteps); //this will vary as velocity will vary with height
printf("Distance advected y = %g\n", vely*dt*(float) nsteps);
/*** Set up Gaussian initial conditions ***/
/* LOOP 3 */
//here we can use collapse keyword as we have a inner loop thats not related to the outerloop by iteration counts
#pragma omp parallel deafult(none) collapse(2) private(x2,y2) //we have to make x and y2 private to not cause unexpected results due to race conditions
for (int i=0; i<NX+2; i++){
for (int j=0; j<NY+2; j++){
x2 = (x-x0) * (x-x0);
y2 = (y[j]-y0) * (y[j]-y0);
u[j] = exp(-1.0 * ( (x2/(2.0*sigmax2)) + (y2/(2.0*sigmay2))));
}
}
/*** Write array of initial u values out to file ***/
FILE *initialfile;
initialfile = fopen("initial.dat", "w");
/* LOOP 4 */
for (int i=0; i<NX+2; i++){ //cannot paralise because this is data output so it has to be delivered in a specific order an multithreading doesnt gurantee that
for (int j=0; j<NY+2; j++){
fprintf(initialfile, "%g %g %g\n", x, y[j], u[j]);
}
}
fclose(initialfile);
/*** Update solution by looping over time steps ***/
/* LOOP 5 */
for (int m=0; m<nsteps; m++){//cant just use collapse here as some loops inside this loop cannot be paralised
#pragma omp parallel deafult(none)
{
/*** Apply boundary conditions at u[0][:] and u[NX+1][:] ***/
/* LOOP 6 */
for (int j=0; j<NY+2; j++){
u[0][j] = bval_left;
u[NX+1][j] = bval_right;
}
/*** Apply boundary conditions at u[:][0] and u[:][NY+1] ***/
/* LOOP 7 */
for (int i=0; i<NX+2; i++){
u[0] = bval_lower;
u[NY+1] = bval_upper;
}
}
/*** Calculate rate of change of u using leftward difference ***/
/* Loop over points in the domain but not boundary values */
/* LOOP 8 */
//cant parallelise because we have a dependency where we need i-1 elements that could of been executed by other threads causing unexpected results
for (int i=1; i<NX+1; i++){
for (int j=1; j<NY+1; j++){
dudt[j] = -velx[j] * (u[j] - u[i-1][j]) / dx //velx is now variable depending on the height
- vely * (u[j] - u[j-1]) / dy;
}
}
/*** Update u from t to t+dt ***/
/* Loop over points in the domain but not boundary values */
/* LOOP 9 */
#pragma omp parallel for deafult(none) collapse(2) reduction(+ : u)//we have to perform a reduction as it would give innacuare results because here we are accumulating a sum so it has to be done in a specific manner for the result to be valid
for (int i=1; i<NX+1; i++){
for (int j=1; j<NY+1; j++){
u[j] = u[j] + dudt[j] * dt[j];
}
}
} // time loop
/*** Write array of final u values out to file ***/
FILE *finalfile;
finalfile = fopen("final.dat", "w");
/* LOOP 10 */
//cant paralise once again because we need ot deliver data in only one order which isnt guranteed by multithreading
for (int i=0; i<NX+2; i++){
for (int j=0; j<NY+2; j++){
fprintf(finalfile, "%g %g %g\n", x, y[j], u[j]);
}
}
fclose(finalfile);
//creating vertical average we want to plot a average value of u for each x point.
FILE *verticalaverage;
verticalaverage = fopen("VerticalAverage.dat", "w");
/* LOOP 11 */
//cant paralise once again because we need ot deliver data in only one order which isnt guranteed by multithreading
float vertical_avg[NX+2][NX+2][NX+2]; //a array that holds average values of u for every x value
int itemCounter = 0 ;
bool found = false;
// {X VALUE , total U VALUE,total number of u values }
//CONSIDR BOUNDARY THO
for (int i=0; i<NX+2; i++){
for (int j=0; j<NY+2; j++){
//IF FOUND ADD U AND AVG
found = false;
for (int k=0; k<itemCounter; k++){ //check if paricular x value was inserted already
if (vertical_avg[k][0][0]==x){
//add to the array a new value of x
vertical_avg[k][1][0]=vertical_avg[k][1][0] + u[j];
vertical_avg[k][1][1]=vertical_avg[k][1][1] + 1; //increments the count as we added a new value
found =true;
break;//exit loop as its done for this iteration
}
}
if (found == false){
//add new x value into the array
vertical_avg[itemCounter][0][0] = x;
vertical_avg[itemCounter][1][0] = u[j];
vertical_avg[itemCounter][1][1] +=1 ;
itemCounter+=1;
}
}
}
//now we need to divide the total sum by number of elements that created the total value of u to find the average u for each x value.
for (int k= 0 ;k<itermcounter;k++){
vertical_avg[itemCounter][1][0] = vertical_avg[itemCounter][1][0]/ vertical_avg[itemCounter][1][1];
//each field here i now the average value of u for every x value.
fprintf(fverticalaverage, "%g %g %g\n", vertical_avg[itemCounter][0][0], vertical_avg[itemCounter][1][0]);
prints x value on the x axis and the u value of the y axis.
}
fclose(verticalaverage);
return 0;
}
