////////////////////////////////////////////////////////////////////////
// THIS FILE IS THE INTERFACE FILE BETWEEN MATLAB AND C //
////////////////////////////////////////////////////////////////////////
// Authors: Matt Ueckermann, Themis Sapsis,
// Pierre Lermusiaux, Pat Haley for MIT course 2.29
////////////////////////////////////
// DECLARE STANDARD HEADER FILES //
////////////////////////////////////
#include "math.h"
#include "mex.h"
#include <stdlib.h>
#include <string.h>
int sign(float v)
{
return((v<0)-(v>0));
}
////////////////////////////
// MASTER FUNCTION FILE //
////////////////////////////
//////////////////////////////////////////
// MATLAB INTERFACE FUNCTION DEFINITION //
//////////////////////////////////////////
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
// Declare input variables
double dx, dy, *ui, *vi, *uj, *vj;
int Nx, Ny, *Nodeu, *Nodev, Nbcs, NbcsDu, NbcsDv;
// Declare output variables
double *Fu, *Fv;
//Declare working variables
int i, j;
double fw, fe, fn, fs, Uiws, Uien, Ujws, Ujen, Vj, Uj;
// Some minor input/output error checking
if (!(nrhs==13)) {
mexErrMsgTxt("Need 13 input arguments");
}
if (!(nlhs==2)) {
mexErrMsgTxt("Need 2 outputs arguments");
}
// Get access to input data from MATLAB
// Pointers are used for large array
// Scalar are re-created in memory
/********!!!! WARNING !!!************/
// WARNING !!!! note that C is expecting to be pointing to type int32 from
// Matlab. Therefore, for this code to work, the C argument must be
// int32. This is done in matlab by int32(C);
// The function is then called like [Bnew Cnew] = axB(5,B,int32(C));
// Use double(Cnew) if you want to convert it back (this might slow things).
/********!!!! WARNING !!!************/
//[Fu Fv] = UVadvect_DO_CDS(Nx, Ny, dx, dy, ui, vi, uj, vj, Nodeu, Nodev, idsu, idsv)
Nx = (int) mxGetScalar( prhs[0] );
Ny = (int) mxGetScalar( prhs[1] );
dx = (double) mxGetScalar( prhs[2] );
dy = (double) mxGetScalar( prhs[3] );
ui = (double*) mxGetPr( prhs[4] );
vi = (double*) mxGetPr( prhs[5] );
uj = (double*) mxGetPr( prhs[6] );
vj = (double*) mxGetPr( prhs[7] );
Nodeu = (int*) mxGetPr( prhs[8] );
Nodev = (int*) mxGetPr( prhs[9] );
Nbcs = (int) mxGetScalar( prhs[10] );
NbcsDu = (int) mxGetScalar( prhs[11] );
NbcsDv = (int) mxGetScalar( prhs[12] );
// idsu = (int*) mxGetPr( prhs[10] );
// Nidsu = mxGetM ( prhs[10] );
// idsv = (int*) mxGetPr( prhs[11] );
// Nidsv = mxGetM ( prhs[11] );
// Example of outputting text to MATLAB
// if (time ==0)
// mexPrintf("Welcome message for time %f\n",time);
// ALLOCATE AND POINT TO OUTPUTS
plhs[0] = mxCreateDoubleMatrix((Nx-1)* Ny , 1, mxREAL);
plhs[1] = mxCreateDoubleMatrix( Nx *(Ny-1), 1, mxREAL);
Fu = (double*)mxGetPr(plhs[0]);
Fv = (double*)mxGetPr(plhs[1]);
for(i=1;i<Nx;i++) {
for(j=1;j<Ny+1;j++) {
if (Nodeu[j +(i)*(Ny+2)]<=Nbcs) //If this is a boundary node, don't do any work
{
Fu[j-1+(i-1)*Ny] = 0.0;
}
else //do work
{
//U advection, quadratic upwind difference scheme
//mexPrintf("%d, %f\n",Nodeu[j +(i-1)*(Ny+2)]-1,dx );
if(Nodeu[j +(i-1)*(Ny+2)]>Nbcs) {
//West face
Uiws = 3.0/8.0*ui[Nodeu[j + i *(Ny+2)]-1] +\
6.0/8.0*ui[Nodeu[j +(i-1)*(Ny+2)]-1] -\
1.0/8.0*ui[Nodeu[j +(i-2)*(Ny+2)]-1];
Uien = 3.0/8.0*ui[Nodeu[j +(i-1)*(Ny+2)]-1] +\
6.0/8.0*ui[Nodeu[j +(i )*(Ny+2)]-1] -\
1.0/8.0*ui[Nodeu[j +(i+1)*(Ny+2)]-1];
Ujws = 3.0/8.0*uj[Nodeu[j + i *(Ny+2)]-1] +\
6.0/8.0*uj[Nodeu[j +(i-1)*(Ny+2)]-1] -\
1.0/8.0*uj[Nodeu[j +(i-2)*(Ny+2)]-1];
Ujen = 3.0/8.0*uj[Nodeu[j +(i-1)*(Ny+2)]-1] +\
6.0/8.0*uj[Nodeu[j +(i )*(Ny+2)]-1] -\
1.0/8.0*uj[Nodeu[j +(i+1)*(Ny+2)]-1];
}
else //treat boundary condition
{
if (Nodeu[j +(i-1)*(Ny+2)]<=NbcsDu) {
//West face
Uiws = 0.5*(ui[Nodeu[j +(i-1)*(Ny+2)]-1]+ui[Nodeu[j +(i )*(Ny+2)]-1]);
Uien = Uiws;
Ujws = 0.5*(uj[Nodeu[j +(i-1)*(Ny+2)]-1]+uj[Nodeu[j +(i )*(Ny+2)]-1]);
Ujen = Ujws;
}
else {
Uiws= (0.5*ui[Nodeu[j +(i-1)*(Ny+2)]-1]+ui[Nodeu[j +(i)*(Ny+2)]-1]);
Uien = Uiws;
Ujws= (0.5*uj[Nodeu[j +(i-1)*(Ny+2)]-1]+uj[Nodeu[j +(i)*(Ny+2)]-1]);
Ujen = Ujws;
}
}
//Treat additional boundary problems. The QUICk scheme won't
// work if neumann boundaries are at i-2 or i+1
if (i-2>=0) { //This is to avoid segfaults
if(Nodeu[j+(i-2)*(Ny+2)]<=Nbcs && Nodeu[j+(i-2)*(Ny+2)]>NbcsDu) {
//if the most west boundary is Neumann use CDS
Uiws = 0.5*(ui[Nodeu[j +(i-1)*(Ny+2)]-1]+ui[Nodeu[j +(i )*(Ny+2)]-1]);
Ujws = 0.5*(uj[Nodeu[j +(i-1)*(Ny+2)]-1]+uj[Nodeu[j +(i )*(Ny+2)]-1]);
}
}
if(Nodeu[j+(i+1)*(Ny+2)]<=Nbcs && Nodeu[j+(i+1)*(Ny+2)]>NbcsDu) {
//if the most east boundary is Neumann use CDS
Uien = 0.5*(ui[Nodeu[j +(i-1)*(Ny+2)]-1]+ui[Nodeu[j +(i )*(Ny+2)]-1]);
Ujen = 0.5*(uj[Nodeu[j +(i-1)*(Ny+2)]-1]+uj[Nodeu[j +(i )*(Ny+2)]-1]);
}
//Do some minor slope limiting
if (i-2>=0) {
if (sign(uj[Nodeu[j +(i-1)*(Ny+2)]-1]) != sign(uj[Nodeu[j +(i-2)*(Ny+2)]-1]))
{
Uiws = ui[Nodeu[j +(i-1)*(Ny+2)]-1];
Ujws = uj[Nodeu[j +(i-1)*(Ny+2)]-1];
}
}
if (sign(uj[Nodeu[j +(i+1)*(Ny+2)]-1]) != sign(uj[Nodeu[j +(i)*(Ny+2)]-1])) {
Uien = ui[Nodeu[j +(i)*(Ny+2)]-1];
Ujen = uj[Nodeu[j +(i)*(Ny+2)]-1];
}
fw = (0.5/dx)*(Uiws*(Ujws+fabs(Ujws)) + Uien*(Ujen-fabs(Ujen)) );
if (Nodeu[j +(i+1)*(Ny+2)]>Nbcs) {
//East face
Uiws = 3.0/8.0*ui[Nodeu[j +(i+1)*(Ny+2)]-1] +\
6.0/8.0*ui[Nodeu[j +(i )*(Ny+2)]-1] -\
1.0/8.0*ui[Nodeu[j +(i-1)*(Ny+2)]-1];
Uien = 3.0/8.0*ui[Nodeu[j +(i )*(Ny+2)]-1] +\
6.0/8.0*ui[Nodeu[j +(i+1)*(Ny+2)]-1] -\
1.0/8.0*ui[Nodeu[j +(i+2)*(Ny+2)]-1];
Ujws = 3.0/8.0*uj[Nodeu[j +(i+1)*(Ny+2)]-1] +\
6.0/8.0*uj[Nodeu[j +(i )*(Ny+2)]-1] -\
1.0/8.0*uj[Nodeu[j +(i-1)*(Ny+2)]-1];
Ujen = 3.0/8.0*uj[Nodeu[j +(i )*(Ny+2)]-1] +\
6.0/8.0*uj[Nodeu[j +(i+1)*(Ny+2)]-1] -\
1.0/8.0*uj[Nodeu[j +(i+2)*(Ny+2)]-1];
}
else //treat boundary condition
{
if (Nodeu[j +(i+1)*(Ny+2)]<=NbcsDu) {
//mexPrintf("Dirichlet %d,%d, %f\n",Nodeu[j+(i+1)*(Ny+2)]-1,NbcsDu,uj[Nodeu[j +(i+1)*(Ny+2)]-1] );
//East face
Uiws = 0.5*(ui[Nodeu[j +(i+1)*(Ny+2)]-1] + ui[Nodeu[j +i*(Ny+2)]-1]);
Uien = Uiws;
Ujws = 0.5*(uj[Nodeu[j +(i+1)*(Ny+2)]-1] + uj[Nodeu[j +i*(Ny+2)]-1]);
Ujen = Ujws;
}
else {
Uiws = (0.5*ui[Nodeu[j +(i+1)*(Ny+2)]-1] + ui[Nodeu[j +i*(Ny+2)]-1]);
Uien = Uiws;
Ujws = (0.5*uj[Nodeu[j +(i+1)*(Ny+2)]-1] + uj[Nodeu[j +i*(Ny+2)]-1]);
Ujen = Ujws;
//mexPrintf("Neumann %d,%d, %f\n",Nodeu[j+(i+1)*(Ny+2)]-1,NbcsDu,uj[Nodeu[j +(i+1)*(Ny+2)]-1] );
}
}
//Treat additional boundary problems. The QUICk scheme won't
// work if neumann boundaries are at i-2 or i+1
if (i + 2 <= Nx) { //This is to avoid segfaults
if(Nodeu[j+(i+2)*(Ny+2)]<=Nbcs && Nodeu[j+(i+2)*(Ny+2)]>NbcsDu) {
//if the most east boundary is Neumann use CDS
Uien = 0.5*(ui[Nodeu[j +(i+1)*(Ny+2)]-1] + ui[Nodeu[j +i*(Ny+2)]-1]);
Ujen = 0.5*(uj[Nodeu[j +(i+1)*(Ny+2)]-1] + uj[Nodeu[j +i*(Ny+2)]-1]);
}
}
if(Nodeu[j+(i-1)*(Ny+2)]<=Nbcs && Nodeu[j+(i-1)*(Ny+2)]>NbcsDu) {
//if the most west boundary is Neumann use CDS
Uiws = 0.5*(ui[Nodeu[j +(i+1)*(Ny+2)]-1] + ui[Nodeu[j +i*(Ny+2)]-1]);
Ujws = 0.5*(uj[Nodeu[j +(i+1)*(Ny+2)]-1] + uj[Nodeu[j +i*(Ny+2)]-1]);
}
//Do some minor slope limiting
if (sign(uj[Nodeu[j +(i-1)*(Ny+2)]-1]) != sign(uj[Nodeu[j +(i)*(Ny+2)]-1])) {
Uiws = ui[Nodeu[j +(i)*(Ny+2)]-1];
Ujws = uj[Nodeu[j +(i)*(Ny+2)]-1];
}
if (i + 2 <= Nx) {
if (sign(uj[Nodeu[j +(i+1)*(Ny+2)]-1]) != sign(uj[Nodeu[j+(i+2)*(Ny+2)]-1])) {
Uien = ui[Nodeu[j +(i + 1)*(Ny+2)]-1];
Ujen = uj[Nodeu[j +(i + 1)*(Ny+2)]-1];
}
}
fe = (0.5/dx)*(Uiws*(Ujws+fabs(Ujws)) + Uien*(Ujen-fabs(Ujen)) );
if (Nodeu[j+1+ i*(Ny+2)]>Nbcs) {
//South Face
Uiws = 3.0/8.0*ui[Nodeu[j + i*(Ny+2)]-1] +\
6.0/8.0*ui[Nodeu[j+1+ i*(Ny+2)]-1] -\
1.0/8.0*ui[Nodeu[j+2+ i*(Ny+2)]-1];
Uien = 3.0/8.0*ui[Nodeu[j+1+ i*(Ny+2)]-1] +\
6.0/8.0*ui[Nodeu[j + i*(Ny+2)]-1] -\
1.0/8.0*ui[Nodeu[j-1+ i*(Ny+2)]-1];
Vj = 0.5*( vj[Nodev[j +(i+1)*(Ny+1)]-1] + vj[Nodev[j +i*(Ny+1)]-1] );
}
else { //treat boundaries
if (Nodeu[j+1+(i)*(Ny+2)]<=NbcsDu) { //Dirichlet on u
Uiws = 0.5*(ui[Nodeu[j+1+ i*(Ny+2)]-1]+ui[Nodeu[j + i*(Ny+2)]-1]);
Uien = Uiws;
}
else { //Neumann on u
Uiws = (0.5*ui[Nodeu[j+1+ i*(Ny+2)]-1]+ui[Nodeu[j + i*(Ny+2)]-1]);
Uien = Uiws;
}
if (Nodev[j +(i+1)*(Ny+1)]<=NbcsDv || Nodev[j +(i+1)*(Ny+1)]>Nbcs) { //Dirichlet on v
Vj = 0.5*( vj[Nodev[j +(i+1)*(Ny+1)]-1] + vj[Nodev[j +i*(Ny+1)]-1] );
}
else {
Vj = 0.5*( vj[Nodev[j +(i+1)*(Ny+1)]-1] + vj[Nodev[j-1 +(i+1)*(Ny+1)]-1] \
+ vj[Nodev[j + i* (Ny+1)]-1] + vj[Nodev[j-1+ i* (Ny+1)]-1] );
}
}
//Treat additional boundary problems. The QUICk scheme won't
// work if neumann boundaries are at i-2 or i+1
if (j+2<=Ny+1) { //This is to avoid segfaults
if(Nodeu[j+2+i*(Ny+2)]<=Nbcs && Nodeu[j+2+i*(Ny+2)]>NbcsDu) {
//if the most north boundary is Neumann use CDS
Uiws = 0.5*(ui[Nodeu[j+1+ i*(Ny+2)]-1]+ui[Nodeu[j + i*(Ny+2)]-1]);
}
}
if(Nodeu[j-1+i*(Ny+2)]<=Nbcs && Nodeu[j-1+i*(Ny+2)]>NbcsDu) {
//if the most south boundary is Neumann use CDS
Uien = 0.5*(ui[Nodeu[j+1+ i*(Ny+2)]-1]+ui[Nodeu[j + i*(Ny+2)]-1]);
}
fs = (0.5/dy)*( Uiws*(Vj+fabs(Vj)) + Uien*(Vj-fabs(Vj)) );
//North Face
if (Nodeu[j-1+ i*(Ny+2)]>Nbcs) {
Uiws = 3.0/8.0*ui[Nodeu[j-1+ i*(Ny+2)]-1] +\
6.0/8.0*ui[Nodeu[j + i*(Ny+2)]-1] -\
1.0/8.0*ui[Nodeu[j+1+ i*(Ny+2)]-1];
Uien = 3.0/8.0*ui[Nodeu[j + i*(Ny+2)]-1] +\
6.0/8.0*ui[Nodeu[j-1+ i*(Ny+2)]-1] -\
1.0/8.0*ui[Nodeu[j-2+ i*(Ny+2)]-1];
Vj = 0.5*( vj[Nodev[j-1 +(i+1)*(Ny+1)]-1] + vj[Nodev[j-1+i*(Ny+1)]-1] );
}
else { //treat boundaries
if (Nodeu[j-1+(i)*(Ny+2)]<=NbcsDu) { //Dirichlet on u
Uiws = 0.5*(ui[Nodeu[j + i*(Ny+2)]-1]+ui[Nodeu[j-1+ i*(Ny+2)]-1]);
Uien = Uiws;
}
else { //Neuman on u
Uiws = (0.5*ui[Nodeu[j + i*(Ny+2)]-1]+ui[Nodeu[j-1+ i*(Ny+2)]-1]);
Uien = Uiws;
}
if (Nodev[j-1 +(i+1)*(Ny+1)]<=NbcsDv || Nodev[j-1 +(i+1)*(Ny+1)]>Nbcs) { //Dirichlet on v
Vj = 0.5*( vj[Nodev[j-1 +(i+1)*(Ny+1)]-1] + vj[Nodev[j-1+i*(Ny+1)]-1] );
}
else { //Neuman on v
Vj = 0.5*( vj[Nodev[j +(i+1)*(Ny+1)]-1] + vj[Nodev[j-1 +(i+1)*(Ny+1)]-1] \
+ vj[Nodev[j + i* (Ny+1)]-1] + vj[Nodev[j-1 + i* (Ny+1)]-1] );
}
}
//Treat additional boundary problems. The QUICk scheme won't
// work if neumann boundaries are at i-2 or i+1
if (j-2>=0) { //This is to avoid segfaults
if(Nodeu[j-2+i*(Ny+2)]<=Nbcs && Nodeu[j-2+i*(Ny+2)]>NbcsDu) {
//if the most north boundary is Neumann use CDS
Uien = 0.5*(ui[Nodeu[j + i*(Ny+2)]-1]+ui[Nodeu[j-1+ i*(Ny+2)]-1]);
}
}
if(Nodeu[j+1+i*(Ny+2)]<=Nbcs && Nodeu[j+1+i*(Ny+2)]>NbcsDu) {
//if the most south boundary is Neumann use CDS
Uiws = 0.5*(ui[Nodeu[j + i*(Ny+2)]-1]+ui[Nodeu[j-1+ i*(Ny+2)]-1]);
}
fn = (0.5/dy)*( Uiws*(Vj+fabs(Vj)) + Uien*(Vj-fabs(Vj)) );
//Calculate Total flux change for this cell
Fu[j-1+(i-1)*Ny] = fw-fe+fs-fn;
}
}
}
for(i=1;i<Nx+1;i++) {
for(j=1;j<Ny;j++) {
if (Nodev[j +(i)*(Ny+1)]<=Nbcs) //If boundary node don't do work
{
Fv[j-1+(i-1)*(Ny-1)] = 0.0;
}
else //do work
{
//V advection, QUICK scheme
//mexPrintf("%d, %f\n",Nodeu[j +(i-1)*(Ny+2)]-1,dx );
if(Nodev[j +(i-1)*(Ny+1)]>Nbcs) {
//West face
Uiws = 3.0/8.0*vi[Nodev[j + i *(Ny+1)]-1] +\
6.0/8.0*vi[Nodev[j +(i-1)*(Ny+1)]-1] -\
1.0/8.0*vi[Nodev[j +(i-2)*(Ny+1)]-1];
Uien = 3.0/8.0*vi[Nodev[j +(i-1)*(Ny+1)]-1] +\
6.0/8.0*vi[Nodev[j +(i )*(Ny+1)]-1] -\
1.0/8.0*vi[Nodev[j +(i+1)*(Ny+1)]-1];
Uj = 0.5*( uj[Nodeu[j+1 +(i-1)*(Ny+2)]-1] + uj[Nodeu[j +(i-1)*(Ny+2)]-1] );
}
else //treat boundary condition
{
//West face
if (Nodev[j +(i-1)*(Ny+1)]-1<=NbcsDv) { //dirichlet on v
Uiws = 0.5*(vi[Nodev[j +(i-1)*(Ny+1)]-1]+vi[Nodev[j +(i )*(Ny+1)]-1]);
Uien = Uiws;
}
else { //Neuman on v
Uiws = (0.5*vi[Nodev[j +(i-1)*(Ny+1)]-1]+vi[Nodev[j +(i )*(Ny+1)]-1]);
Uien = Uiws;
}
if (Nodeu[j+1 +(i-1)*(Ny+2)]<=NbcsDu || Nodeu[j+1 +(i-1)*(Ny+2)]>Nbcs) { //dirichlet on u
Uj = 0.5*( uj[Nodeu[j+1 +(i-1)*(Ny+2)]-1] + uj[Nodeu[j +(i-1)*(Ny+2)]-1] );
}
else {
Uj = 0.5*( uj[Nodeu[j+1 +(i-1)*(Ny+2)]-1]+uj[Nodeu[j+1 +(i)*(Ny+2)]-1]\
+ uj[Nodeu[j + (i-1)*(Ny+2)]-1]+uj[Nodeu[j +(i)*(Ny+2)]-1] );
}
}
//Treat additional boundary problems. The QUICk scheme won't
// work if neumann boundaries are at i-2 or i+1
if (i-2>=0) { //This is to avoid segfaults
if(Nodev[j+(i-2)*(Ny+1)]<=Nbcs && Nodev[j+(i-2)*(Ny+1)]>NbcsDv) {
//if the most west boundary is Neumann use CDS
Uiws = 0.5*(vi[Nodev[j +(i-1)*(Ny+1)]-1]+vi[Nodev[j +(i )*(Ny+1)]-1]);
}
}
if(Nodev[j+(i+1)*(Ny+1)]<=Nbcs && Nodev[j+(i+1)*(Ny+1)]>NbcsDv) {
//if the most east boundary is Neumann use CDS
Uien = 0.5*(vi[Nodev[j +(i-1)*(Ny+1)]-1]+vi[Nodev[j +(i )*(Ny+1)]-1]);
}
fw = (0.5/dx)*(Uiws*(Uj+fabs(Uj)) + Uien*(Uj-fabs(Uj)) );
if (Nodev[j +(i+1)*(Ny+1)]>Nbcs) {
//East face
Uiws = 3.0/8.0*vi[Nodev[j +(i+1)*(Ny+1)]-1] +\
6.0/8.0*vi[Nodev[j +(i )*(Ny+1)]-1] -\
1.0/8.0*vi[Nodev[j +(i-1)*(Ny+1)]-1];
Uien = 3.0/8.0*vi[Nodev[j +(i )*(Ny+1)]-1] +\
6.0/8.0*vi[Nodev[j +(i+1)*(Ny+1)]-1] -\
1.0/8.0*vi[Nodev[j +(i+2)*(Ny+1)]-1];
Uj = 0.5*( uj[Nodeu[j+1 + i*(Ny+2) ]-1] + uj[Nodeu[j +i*(Ny+2)]-1] );
}
else //treat boundary condition
{
//East face
if (Nodev[j +(i+1)*(Ny+1)]-1<=NbcsDv) {//Dirichlet on v
Uiws = 0.5*(vi[Nodev[j +(i )*(Ny+1)]-1]+vi[Nodev[j +(i+1)*(Ny+1)]-1]);
Uien = Uiws;
}
else { //Neumann on v
Uiws = (0.5*vi[Nodev[j +(i )*(Ny+1)]-1]+vi[Nodev[j +(i+1)*(Ny+1)]-1]);
Uien = Uiws;
}
if (Nodeu[j+1 +(i)*(Ny+2)]<=NbcsDu || Nodeu[j+1 +(i)*(Ny+2)]>Nbcs) { //dirichlet on u
Uj = 0.5*( uj[Nodeu[j+1 + i*(Ny+2) ]-1] + uj[Nodeu[j +i*(Ny+2)]-1] );
}
else { //neuman on u
Uj = 0.5*( uj[Nodeu[j+1 +(i-1)*(Ny+2)]-1]+uj[Nodeu[j+1 +(i)*(Ny+2)]-1]\
+ uj[Nodeu[j + (i-1)*(Ny+2)]-1]+uj[Nodeu[j +(i)*(Ny+2)]-1] );
}
}
//Treat additional boundary problems. The QUICk scheme won't
// work if neumann boundaries are at i-2 or i+1
if (i+2<=Nx+1) { //This is to avoid segfaults
if(Nodev[j+(i+2)*(Ny+1)]<=Nbcs && Nodev[j+(i+2)*(Ny+1)]>NbcsDv) {
//if the most east boundary is Neumann use CDS
Uien = 0.5*(vi[Nodev[j +(i )*(Ny+1)]-1]+vi[Nodev[j +(i+1)*(Ny+1)]-1]);
}
}
if(Nodev[j+(i-1)*(Ny+1)]<=Nbcs && Nodev[j+(i-1)*(Ny+1)]>NbcsDv) {
//if the most west boundary is Neumann use CDS
Uiws = 0.5*(vi[Nodev[j +(i )*(Ny+1)]-1]+vi[Nodev[j +(i+1)*(Ny+1)]-1]);
}
fe = (0.5/dx)*(Uiws*(Uj+fabs(Uj)) + Uien*(Uj-fabs(Uj)) );
if (Nodev[j+1+ i*(Ny+1)]>Nbcs) {
//South Face
Uiws = 3.0/8.0*vi[Nodev[j + i*(Ny+1)]-1] +\
6.0/8.0*vi[Nodev[j+1+ i*(Ny+1)]-1] -\
1.0/8.0*vi[Nodev[j+2+ i*(Ny+1)]-1];
Uien = 3.0/8.0*vi[Nodev[j+1+ i*(Ny+1)]-1] +\
6.0/8.0*vi[Nodev[j + i*(Ny+1)]-1] -\
1.0/8.0*vi[Nodev[j-1+ i*(Ny+1)]-1];
Ujws = 3.0/8.0*vj[Nodev[j + i*(Ny+1)]-1] +\
6.0/8.0*vj[Nodev[j+1+ i*(Ny+1)]-1] -\
1.0/8.0*vj[Nodev[j+2+ i*(Ny+1)]-1];
Ujen = 3.0/8.0*vj[Nodev[j+1+ i*(Ny+1)]-1] +\
6.0/8.0*vj[Nodev[j + i*(Ny+1)]-1] -\
1.0/8.0*vj[Nodev[j-1+ i*(Ny+1)]-1];
}
else {
if (Nodev[j+1+(i)*(Ny+1)]-1<=NbcsDv) {
Uiws = 0.5*(vi[Nodev[j+1+ i*(Ny+1)]-1]+vi[Nodev[j + i*(Ny+1)]-1]);
Uien = Uiws;
Ujws = 0.5*(vj[Nodev[j+1+ i*(Ny+1)]-1]+vj[Nodev[j + i*(Ny+1)]-1]);
Ujen = Ujws;
}
else {
Uiws = (0.5*vi[Nodev[j+1+ i*(Ny+1)]-1]+vi[Nodev[j + i*(Ny+1)]-1]);
Uien = Uiws;
Ujws = (0.5*vj[Nodev[j+1+ i*(Ny+1)]-1]+vj[Nodev[j + i*(Ny+1)]-1]);
Ujen = Ujws;
}
}
//Treat additional boundary problems. The QUICk scheme won't
// work if neumann boundaries are at i-2 or i+1
if (j+2<=Ny) { //This is to avoid segfaults
if(Nodev[j+2+i*(Ny+1)]<=Nbcs && Nodev[j+2+i*(Ny+1)]>NbcsDv) {
//if the most south boundary is Neumann use CDS
Uiws = 0.5*(vi[Nodev[j+1+ i*(Ny+1)]-1]+vi[Nodev[j + i*(Ny+1)]-1]);
Ujws = 0.5*(vj[Nodev[j+1+ i*(Ny+1)]-1]+vj[Nodev[j + i*(Ny+1)]-1]);
}
}
if(Nodev[j-1+i*(Ny+1)]<=Nbcs && Nodev[j-1+i*(Ny+1)]>NbcsDv) {
//if the most north boundary is Neumann use CDS
Uien = 0.5*(vi[Nodev[j+1+ i*(Ny+1)]-1]+vi[Nodev[j + i*(Ny+1)]-1]);
Ujen = 0.5*(vj[Nodev[j+1+ i*(Ny+1)]-1]+vj[Nodev[j + i*(Ny+1)]-1]);
}
//Do some minor slope limiting
if (j+2<=Ny) {
if (sign(vj[Nodev[j+1 + i*(Ny+1)]-1]) != sign(vj[Nodev[j+2 + i*(Ny+1)]-1]))
{
Uiws = vi[Nodev[j+1 +(i)*(Ny+1)]-1];
Ujws = vj[Nodev[j+1 +(i)*(Ny+1)]-1];
}
}
if (sign(vj[Nodev[j-1 + i*(Ny+1)]-1]) != sign(vj[Nodev[j + i*(Ny+1)]-1])) {
Uien = vi[Nodev[j + i*(Ny+1)]-1];
Ujen = vj[Nodev[j + i*(Ny+1)]-1];
}
fs = (0.5/dy)*( Uiws*(Ujws+fabs(Ujws)) + Uien*(Ujen-fabs(Ujen)) );
//North Face
if (Nodev[j-1+ i*(Ny+1)]>Nbcs) {
Uiws = 3.0/8.0*vi[Nodev[j-1+ i*(Ny+1)]-1] +\
6.0/8.0*vi[Nodev[j + i*(Ny+1)]-1] -\
1.0/8.0*vi[Nodev[j+1+ i*(Ny+1)]-1];
Uien = 3.0/8.0*vi[Nodev[j + i*(Ny+1)]-1] +\
6.0/8.0*vi[Nodev[j-1+ i*(Ny+1)]-1] -\
1.0/8.0*vi[Nodev[j-2+ i*(Ny+1)]-1];
Ujws = 3.0/8.0*vj[Nodev[j-1+ i*(Ny+1)]-1] +\
6.0/8.0*vj[Nodev[j + i*(Ny+1)]-1] -\
1.0/8.0*vj[Nodev[j+1+ i*(Ny+1)]-1];
Ujen = 3.0/8.0*vj[Nodev[j + i*(Ny+1)]-1] +\
6.0/8.0*vj[Nodev[j-1+ i*(Ny+1)]-1] -\
1.0/8.0*vj[Nodev[j-2+ i*(Ny+1)]-1];
}
else {
if (Nodev[j-1+(i)*(Ny+1)]-1<=NbcsDv) {
Uiws = 0.5*(vi[Nodev[j + i*(Ny+1)]-1]+vi[Nodev[j-1+ i*(Ny+1)]-1]);
Uien = Uiws;
Ujws = 0.5*(vj[Nodev[j + i*(Ny+1)]-1]+vj[Nodev[j-1+ i*(Ny+1)]-1]);
Ujen = Ujws;
}
else {
Uiws = (0.5*vi[Nodev[j + i*(Ny+1)]-1]+vi[Nodev[j-1+ i*(Ny+1)]-1]);
Uien = Uiws;
Ujws = (0.5*vj[Nodev[j + i*(Ny+1)]-1]+vj[Nodev[j-1+ i*(Ny+1)]-1]);
Ujen = Ujws;
}
}
//Treat additional boundary problems. The QUICk scheme won't
// work if neumann boundaries are at i-2 or i+1
if (j-2>=0) { //This is to avoid segfaults
if(Nodev[j-2+i*(Ny+1)]<=Nbcs && Nodev[j-2+i*(Ny+1)]>NbcsDv) {
//if the most north boundary is Neumann use CDS
Uien = 0.5*(vi[Nodev[j + i*(Ny+1)]-1]+vi[Nodev[j-1+ i*(Ny+1)]-1]);
Ujen = 0.5*(vj[Nodev[j + i*(Ny+1)]-1]+vj[Nodev[j-1+ i*(Ny+1)]-1]);
}
}
if(Nodev[j+1+i*(Ny+1)]<=Nbcs && Nodev[j+1+i*(Ny+1)]>NbcsDv) {
//if the most south boundary is Neumann use CDS
Uiws = 0.5*(vi[Nodev[j + i*(Ny+1)]-1]+vi[Nodev[j-1+ i*(Ny+1)]-1]);
Ujws = 0.5*(vj[Nodev[j + i*(Ny+1)]-1]+vj[Nodev[j-1+ i*(Ny+1)]-1]);
}
//Do some minor slope limiting
if (sign(vj[Nodev[j+1 + i*(Ny+1)]-1]) != sign(vj[Nodev[j + i*(Ny+1)]-1])) {
Uiws = vi[Nodev[j + i*(Ny+1)]-1];
Ujws = vj[Nodev[j + i*(Ny+1)]-1];
}
if (j-2>=0) {
if (sign(vj[Nodev[j-1 + i*(Ny+1)]-1]) != sign(vj[Nodev[j-2 + i*(Ny+1)]-1])) {
Uien = vi[Nodev[j-1 + i*(Ny+1)]-1];
Ujen = vj[Nodev[j-1 + i*(Ny+1)]-1];
}
}
fn = (0.5/dy)*( Uiws*(Ujws+fabs(Ujws)) + Uien*(Ujen-fabs(Ujen)) );
Fv[j-1+(i-1)*(Ny-1)] = fw-fe+fs-fn;
}
}
}
}