////////////////////////////////////////////////////////////////////////
// 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>
#define MIN(a,b) ((a)>(b)?(b):(a))
#define MAX(a,b) ((a)<(b)?(b):(a))
int sign(float v)
{
return((v<0)-(v>0));
}
/////////////////////////////
// DEFINE LIMITER FUNCTION //
/////////////////////////////
double TVD(double phi)
{
//MC limiter
return(MAX(0.0, MIN(MIN(0.5 * (1.0 + phi), 2.0), 2.0 * phi)));
//minmod limiter
// return(((phi) > (0) ? MIN(phi,1) : 0.0));
//superbee
// return(MAX(MAX(0,MIN(1.0, 2.0 * phi)),MIN(2,phi)));
//Upwind
// return(0.0);
// Lax Wendroff (NOT TVD)
// return(1.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, dt, *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 fl, fr, fu, fd, theta, Uj, sumQ, diffQ;
// Some minor input/output error checking
if (!(nrhs==14)) {
mexErrMsgTxt("Need 14 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] );
dt = (double) mxGetScalar( prhs[4] );
ui = (double*) mxGetPr( prhs[5] );
vi = (double*) mxGetPr( prhs[6] );
uj = (double*) mxGetPr( prhs[7] );
vj = (double*) mxGetPr( prhs[8] );
Nodeu = (int*) mxGetPr( prhs[9] );
Nodev = (int*) mxGetPr( prhs[10] );
Nbcs = (int) mxGetScalar( prhs[11] );
NbcsDu = (int) mxGetScalar( prhs[12] );
NbcsDv = (int) mxGetScalar( prhs[13] );
// 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);
// These macros allow easy access to values of u and v at
// UU
// U
// LL L X R RR
// D
// DD
// with respect to the current cell
//First define macros for U
#define UM ui[Nodeu[j + i * (Ny + 2)] - 1]
#define UUU ui[Nodeu[j - 2 + i * (Ny + 2)] - 1]
#define UU ui[Nodeu[j - 1 + i *(Ny + 2)] - 1]
#define UD ui[Nodeu[j + 1 + i * (Ny + 2)] - 1]
#define UDD ui[Nodeu[j + 2 + i * (Ny + 2)] - 1]
#define ULL ui[Nodeu[j + (i - 2) * (Ny + 2)] - 1]
#define UL ui[Nodeu[j + (i - 1) * (Ny + 2)] - 1]
#define UR ui[Nodeu[j + (i + 1) * (Ny + 2)] - 1]
#define URR ui[Nodeu[j + (i + 2) * (Ny + 2)] - 1]
//We will also need to access an averaged value of V at the cell
//interface
#define VAVEU 0.5 * (vj[Nodev[j - 1 +(i + 1) * (Ny + 1)] - 1] + vj[Nodev[j - 1 + i * (Ny + 1)] - 1] )
#define VAVED 0.5 * (vj[Nodev[j +(i + 1) * (Ny + 1)] - 1] + vj[Nodev[j + i * (Ny + 1)] - 1] )
#define UjL 0.5 * (uj[Nodeu[j + (i - 1) * (Ny + 2)] - 1] + uj[Nodeu[j + (i) * (Ny + 2)] - 1]) //\
+0.5 * (fabs(uj[Nodeu[j + (i - 1) * (Ny + 2)] - 1]) - fabs(uj[Nodeu[j + (i) * (Ny + 2)] - 1]))
#define UjR 0.5 * (uj[Nodeu[j + (i + 1) * (Ny + 2)] - 1] + uj[Nodeu[j + (i) * (Ny + 2)] - 1]) //\
-0.5 * (fabs(uj[Nodeu[j + (i + 1) * (Ny + 2)] - 1]) - fabs(uj[Nodeu[j + (i) * (Ny + 2)] - 1]))
//Now Define macros for V
#define VM vi[Nodev[j + i * (Ny + 1)] - 1]
#define VUU vi[Nodev[j - 2 + i * (Ny + 1)] - 1]
#define VU vi[Nodev[j - 1 + i *(Ny + 1)] - 1]
#define VD vi[Nodev[j + 1 + i * (Ny + 1)] - 1]
#define VDD vi[Nodev[j + 2 + i * (Ny + 1)] - 1]
#define VLL vi[Nodev[j + (i - 2) * (Ny + 1)] - 1]
#define VL vi[Nodev[j + (i - 1) * (Ny + 1)] - 1]
#define VR vi[Nodev[j + (i + 1) * (Ny + 1)] - 1]
#define VRR vi[Nodev[j + (i + 2) * (Ny + 1)] - 1]
//We will also need to access an averaged value of V at the cell
//interface
#define UAVEL 0.5 * ( uj[Nodeu[j + 1 +(i - 1) * (Ny + 2)] -1 ] + uj[Nodeu[j + (i - 1) * (Ny + 2)] - 1] )
#define UAVER 0.5 * ( uj[Nodeu[j + 1 + i * (Ny + 2) ] - 1] + uj[Nodeu[j + i * (Ny + 2)] - 1] )
#define VjD 0.5 * (vj[Nodev[j + i * (Ny + 1)] - 1] + vj[Nodev[j + 1 + i * (Ny + 1)] - 1]) //\
-0.5 * (fabs(vj[Nodev[j + i * (Ny + 1)] - 1]) - fabs(vj[Nodev[j + 1 + i * (Ny + 1)] - 1]))
#define VjU 0.5 * (vj[Nodev[j + i * (Ny + 1)] - 1] + vj[Nodev[j - 1 + i * (Ny + 1)] - 1]) //\
+0.5 * (fabs(vj[Nodev[j + i * (Ny + 1)] - 1]) - fabs(vj[Nodev[j - 1 + i * (Ny + 1)] - 1]))
// 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
sumQ = UM + UL;
diffQ = UM - UL;
Uj = UjL;
theta = ((Uj + fabs(Uj)) / 2.0 * (UL - ULL) \
+ (Uj - fabs(Uj)) / 2.0 * (UR - UM)) \
/ (diffQ) / Uj;
}
else //treat boundary condition
{
if (Nodeu[j +(i-1)*(Ny+2)]<=NbcsDu) {
//West face
sumQ = UM + UL;
diffQ = UM - UL;
Uj = UjL;
theta = ((Uj + fabs(Uj)) / 2.0 * (diffQ) \
+ (Uj - fabs(Uj)) / 2.0 * (UR - UM)) \
/ (diffQ) / Uj;
}
else {
sumQ = 2.0 * UM + UL;
diffQ = UM - (UM + UL);
Uj = 2.0 * UjL;
theta = ((Uj + fabs(Uj)) / 2.0 * (diffQ) \
+ (Uj - fabs(Uj)) / 2.0 * (UR - UM)) \
/ (diffQ) / Uj;
}
}
//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:
theta = ((Uj + fabs(Uj)) / 2.0 * (UL - (ULL + UL)) \
+ (Uj - fabs(Uj)) / 2.0 * (UR - UM)) \
/ (diffQ) / Uj;
}
}
if(Nodeu[j+(i+1)*(Ny+2)]<=Nbcs && Nodeu[j+(i+1)*(Ny+2)]>NbcsDu) {
//if the most east boundary is Neumann use:
theta = ((Uj + fabs(Uj)) / 2.0 * (UL - ULL) \
+ (Uj - fabs(Uj)) / 2.0 * ((UM + UR) - UM)) \
/ (diffQ) / Uj;
}
fl = 0.5 * Uj * (sumQ) \
+ 0.5 * fabs(Uj) * (-1.0 + (1.0 - fabs(Uj * dt / dx)) * TVD(theta)) \
* (diffQ);
if (Nodeu[j +(i+1)*(Ny+2)]>Nbcs) {
//East face
sumQ = UR + UM;
diffQ = UR - UM;
Uj = UjR;
theta = ((Uj + fabs(Uj)) / 2.0 * (UM - UL) \
+ (Uj - fabs(Uj)) / 2.0 * (URR - UR)) \
/ (diffQ) / Uj;
}
else //treat boundary condition
{
if (Nodeu[j +(i+1)*(Ny+2)]<=NbcsDu) {
//East face
sumQ = UR + UM;
diffQ = UR - UM;
Uj = UjR;
theta = ((Uj + fabs(Uj)) / 2.0 * (UM - UL) \
+ (Uj - fabs(Uj)) / 2.0 * (diffQ)) \
/ (diffQ) / Uj;
}
else {
sumQ = UR + 2.0 * UM;
diffQ = (UM + UR) - UM;
Uj = 2.0 * UjR;
theta = ((Uj + fabs(Uj)) / 2.0 * (UM - UL) \
+ (Uj - fabs(Uj)) / 2.0 * (diffQ)) \
/ (diffQ) / Uj;
}
}
//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
theta = ((Uj + fabs(Uj)) / 2.0 * (UM - UL) \
+ (Uj - fabs(Uj)) / 2.0 * ((UR + URR) - UR)) \
/ (diffQ) / Uj;
}
}
if(Nodeu[j+(i-1)*(Ny+2)]<=Nbcs && Nodeu[j+(i-1)*(Ny+2)]>NbcsDu) {
//if the most west boundary is Neumann use CDS
theta = ((Uj + fabs(Uj)) / 2.0 * (UM - (UM + UL)) \
+ (Uj - fabs(Uj)) / 2.0 * (URR - UR)) \
/ (diffQ) / Uj;
}
fr = 0.5 * Uj * (sumQ) \
+ 0.5 * fabs(Uj) * (-1.0 + (1.0 - fabs(Uj * dt / dx)) * TVD(theta)) \
* (diffQ);
if (Nodeu[j+1+ i*(Ny+2)]>Nbcs) {
//South Face
sumQ = UM + UD;
diffQ = UM - UD;
Uj = VAVED;
theta = ((Uj + fabs(Uj)) / 2.0 * (UD - UDD) \
+ (Uj - fabs(Uj)) / 2.0 * (UU - UM)) \
/ (diffQ) / Uj;
}
else { //treat boundaries
if (Nodev[j +(i+1)*(Ny+1)]<=NbcsDv || Nodev[j +(i+1)*(Ny+1)]>Nbcs) { //Dirichlet on v
Uj = VAVED;
}
else {
Uj = 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] );
}
if (Nodeu[j+1+(i)*(Ny+2)]<=NbcsDu) { //Dirichlet on u
sumQ = UM + UD;
diffQ = UM - UD;
theta = ((Uj + fabs(Uj)) / 2.0 * (diffQ) \
+ (Uj - fabs(Uj)) / 2.0 * (UU - UM)) \
/ (diffQ) / Uj;
}
else { //Neumann on u
sumQ = 2.0 * UM + UD;
diffQ = UM - (UM + UD);
theta = ((Uj + fabs(Uj)) / 2.0 * (diffQ) \
+ (Uj - fabs(Uj)) / 2.0 * (UU - UM)) \
/ (diffQ) / Uj;
}
}
//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
theta = ((Uj + fabs(Uj)) / 2.0 * (UD - (UD + UDD)) \
+ (Uj - fabs(Uj)) / 2.0 * (UU - UM)) \
/ (diffQ) / Uj;
}
}
if(Nodeu[j-1+i*(Ny+2)]<=Nbcs && Nodeu[j-1+i*(Ny+2)]>NbcsDu) {
//if the most south boundary is Neumann use CDS
theta = ((Uj + fabs(Uj)) / 2.0 * (UD - UDD) \
+ (Uj - fabs(Uj)) / 2.0 * ((UU + UM) - UM)) \
/ (diffQ) / Uj;
}
fd = 0.5 * Uj * (sumQ) \
+ 0.5 * fabs(Uj) * (-1.0 + (1.0 - fabs(Uj * dt / dy)) * TVD(theta)) \
* (diffQ);
//North Face
if (Nodeu[j-1+ i*(Ny+2)]>Nbcs) {
sumQ = UU + UM;
diffQ = UU - UM;
Uj = VAVEU;
theta = ((Uj + fabs(Uj)) / 2.0 * (UM - UD) \
+ (Uj - fabs(Uj)) / 2.0 * (UUU - UU)) \
/ (diffQ) / Uj;
}
else { //treat boundaries
if (Nodev[j-1 +(i+1)*(Ny+1)]<=NbcsDv || Nodev[j-1 +(i+1)*(Ny+1)]>Nbcs) { //Dirichlet on v
Uj = VAVEU;
}
else { //Neuman on v
Uj = 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] );
}
if (Nodeu[j-1+(i)*(Ny+2)]<=NbcsDu) { //Dirichlet on u
sumQ = UU + UM;
diffQ = UU - UM;
theta = ((Uj + fabs(Uj)) / 2.0 * (UM - UD) \
+ (Uj - fabs(Uj)) / 2.0 * (diffQ)) \
/ (diffQ) / Uj;
}
else { //Neuman on u
sumQ = 2.0 * UM + UU;
diffQ = UM - (UM + UU);
theta = ((Uj + fabs(Uj)) / 2.0 * (UM - UD) \
+ (Uj - fabs(Uj)) / 2.0 * (diffQ)) \
/ (diffQ) / Uj;
}
}
//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
theta = ((Uj + fabs(Uj)) / 2.0 * (UM - UD) \
+ (Uj - fabs(Uj)) / 2.0 * ((UUU + UU) - UU)) \
/ (diffQ) / Uj;
}
}
if(Nodeu[j+1+i*(Ny+2)]<=Nbcs && Nodeu[j+1+i*(Ny+2)]>NbcsDu) {
//if the most south boundary is Neumann use CDS
theta = ((Uj + fabs(Uj)) / 2.0 * (UM - (UD + UM)) \
+ (Uj - fabs(Uj)) / 2.0 * (UUU - UU)) \
/ (diffQ) / Uj;
}
fu = 0.5 * Uj * (sumQ) \
+ 0.5 * fabs(Uj) * (-1.0 + (1.0 - fabs(Uj * dt / dy)) * TVD(theta)) \
* (diffQ);
//Calculate Total flux change for this cell
Fu[j-1+(i-1)*Ny] = (fl - fr) / dx + (fd - fu) / dy;
}
}
}
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
sumQ = VM + VL;
diffQ = VM - VL;
Uj = UAVEL;
theta = ((Uj + fabs(Uj)) / 2.0 * (VL - VLL) \
+ (Uj - fabs(Uj)) / 2.0 * (VR - VM)) \
/ (diffQ) / Uj;
}
else //treat boundary condition
{
if (Nodeu[j+1 +(i-1)*(Ny+2)]<=NbcsDu || Nodeu[j+1 +(i-1)*(Ny+2)]>Nbcs) { //dirichlet on u
Uj = UAVEL;
}
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] );
}
//West face
if (Nodev[j +(i-1)*(Ny+1)]-1<=NbcsDv) { //dirichlet on v
sumQ = VM + VL;
diffQ = VM - VL;
theta = ((Uj + fabs(Uj)) / 2.0 * (diffQ) \
+ (Uj - fabs(Uj)) / 2.0 * (VR - VM)) \
/ (diffQ) / Uj;
}
else { //Neuman on v
sumQ = 2.0 * VM + VL;
diffQ = VM - (VL + VM);
theta = ((Uj + fabs(Uj)) / 2.0 * (diffQ) \
+ (Uj - fabs(Uj)) / 2.0 * (VR - VM)) \
/ (diffQ) / Uj;
}
}
//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
theta = ((Uj + fabs(Uj)) / 2.0 * (VL - (VL + VLL)) \
+ (Uj - fabs(Uj)) / 2.0 * (VR - VM)) \
/ (diffQ) / Uj;
}
}
if(Nodev[j+(i+1)*(Ny+1)]<=Nbcs && Nodev[j+(i+1)*(Ny+1)]>NbcsDv) {
//if the most east boundary is Neumann use CDS
theta = ((Uj + fabs(Uj)) / 2.0 * (VL - VLL) \
+ (Uj - fabs(Uj)) / 2.0 * ((VR + VM) - VM)) \
/ (diffQ) / Uj;
}
fl = 0.5 * Uj * (sumQ) \
+ 0.5 * fabs(Uj) * (-1.0 + (1.0 - fabs(Uj * dt / dx)) * TVD(theta)) \
* (diffQ);
if (Nodev[j +(i+1)*(Ny+1)]>Nbcs) {
//East face
sumQ = VR + VM;
diffQ = VR - VM;
Uj = UAVER;
theta = ((Uj + fabs(Uj)) / 2.0 * (VM - VL) \
+ (Uj - fabs(Uj)) / 2.0 * (VRR - VR)) \
/ (diffQ) / Uj;
}
else //treat boundary condition
{
//East face
if (Nodeu[j+1 +(i)*(Ny+2)]<=NbcsDu || Nodeu[j+1 +(i)*(Ny+2)]>Nbcs) { //dirichlet on u
Uj = UAVER;
}
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] );
}
if (Nodev[j +(i+1)*(Ny+1)]-1<=NbcsDv) {//Dirichlet on v
sumQ = VR + VM;
diffQ = VR - VM;
theta = ((Uj + fabs(Uj)) / 2.0 * (VM - VL) \
+ (Uj - fabs(Uj)) / 2.0 * (diffQ)) \
/ (diffQ) / Uj;
}
else { //Neumann on v
sumQ = VR + 2.0 * VM;
diffQ = (VM + VR) - VM;
theta = ((Uj + fabs(Uj)) / 2.0 * (VM - VL) \
+ (Uj - fabs(Uj)) / 2.0 * (diffQ)) \
/ (diffQ) / Uj;
}
}
//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
theta = ((Uj + fabs(Uj)) / 2.0 * (VM - VL) \
+ (Uj - fabs(Uj)) / 2.0 * ((VRR + VR) - VR)) \
/ (diffQ) / Uj;
}
}
if(Nodev[j+(i-1)*(Ny+1)]<=Nbcs && Nodev[j+(i-1)*(Ny+1)]>NbcsDv) {
//if the most west boundary is Neumann use CDS
theta = ((Uj + fabs(Uj)) / 2.0 * (VM - (VL + VM)) \
+ (Uj - fabs(Uj)) / 2.0 * (VRR - VR)) \
/ (diffQ) / Uj;
}
fr = 0.5 * Uj * (sumQ) \
+ 0.5 * fabs(Uj) * (-1.0 + (1.0 - fabs(Uj * dt / dx)) * TVD(theta)) \
* (diffQ);
if (Nodev[j+1+ i*(Ny+1)]>Nbcs) {
//South Face
sumQ = VM + VD;
diffQ = VM - VD;
Uj = VjD;
theta = ((Uj + fabs(Uj)) / 2.0 * (VD - VDD) \
+ (Uj - fabs(Uj)) / 2.0 * (VU - VM)) \
/ (diffQ) / Uj;
}
else {
if (Nodev[j+1+(i)*(Ny+1)]-1<=NbcsDv) {
sumQ = VM + VD;
diffQ = VM - VD;
Uj = VjD;
theta = ((Uj + fabs(Uj)) / 2.0 * (diffQ) \
+ (Uj - fabs(Uj)) / 2.0 * (VU - VM)) \
/ (diffQ) / Uj;
}
else {
sumQ = 2.0 * VM + VD;
diffQ = VM - (VM + VD);
Uj = 2.0 * VjD;
theta = ((Uj + fabs(Uj)) / 2.0 * (diffQ) \
+ (Uj - fabs(Uj)) / 2.0 * (VU - VM)) \
/ (diffQ) / Uj;
}
}
//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
theta = ((Uj + fabs(Uj)) / 2.0 * (VD - ( VD + VDD)) \
+ (Uj - fabs(Uj)) / 2.0 * (VU - VM)) \
/ (diffQ) / Uj;
}
}
if(Nodev[j-1+i*(Ny+1)]<=Nbcs && Nodev[j-1+i*(Ny+1)]>NbcsDv) {
//if the most north boundary is Neumann use CDS
theta = ((Uj + fabs(Uj)) / 2.0 * (VD - VDD) \
+ (Uj - fabs(Uj)) / 2.0 * ((VU + VM) - VM)) \
/ (diffQ) / Uj;
}
fd = 0.5 * Uj * (sumQ) \
+ 0.5 * fabs(Uj) * (-1.0 + (1.0 - fabs(Uj * dt / dy)) * TVD(theta)) \
* (diffQ);
//North Face
if (Nodev[j-1+ i*(Ny+1)]>Nbcs) {
sumQ = VU + VM;
diffQ = VU - VM;
Uj = VjU;
theta = ((Uj + fabs(Uj)) / 2.0 * (VM - VD) \
+ (Uj - fabs(Uj)) / 2.0 * (VUU - VU)) \
/ (diffQ) / Uj;
}
else {
if (Nodev[j-1+(i)*(Ny+1)]-1<=NbcsDv) {
sumQ = VU + VM;
diffQ = VU - VM;
Uj = VjU;
theta = ((Uj + fabs(Uj)) / 2.0 * (VM - VD) \
+ (Uj - fabs(Uj)) / 2.0 * (diffQ)) \
/ (diffQ) / Uj;
}
else {
sumQ = VU + 2.0 * VM;
diffQ = (VU + VM) - VM;
Uj = 2.0 * VjU;
theta = ((Uj + fabs(Uj)) / 2.0 * (VM - VD) \
+ (Uj - fabs(Uj)) / 2.0 * (diffQ)) \
/ (diffQ) / Uj;
}
}
//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
theta = ((Uj + fabs(Uj)) / 2.0 * (VM - VD) \
+ (Uj - fabs(Uj)) / 2.0 * ((VUU + VU) - VU)) \
/ (diffQ) / Uj;
}
}
if(Nodev[j+1+i*(Ny+1)]<=Nbcs && Nodev[j+1+i*(Ny+1)]>NbcsDv) {
//if the most south boundary is Neumann use CDS
theta = ((Uj + fabs(Uj)) / 2.0 * (VM - (VM + VD)) \
+ (Uj - fabs(Uj)) / 2.0 * (VUU - VU)) \
/ (diffQ) / Uj;
}
fu = 0.5 * Uj * (sumQ) \
+ 0.5 * fabs(Uj) * (-1.0 + (1.0 - fabs(Uj * dt / dy)) * TVD(theta)) \
* (diffQ);
//Calculate total flux for this cell
Fv[j-1+(i-1)*(Ny-1)] = (fl - fr) / dx + (fd - fu) / dy;
}
}
}
}