////////////////////////////////////////////////////////////////////////
// THIS FILE IS THE INTERFACE FILE BETWEEN MATLAB AND C //
////////////////////////////////////////////////////////////////////////
// Authors: Matt Ueckermann, 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, NBlen, *u, *v, *rho;
int Nx, Ny, *Nodeu, *Nodev, *NodeP, Nbcs, NbcsDP, NbcsDu, NbcsDv;
// Declare output variables
double *Fu;
//Declare working variables
int i, j;
double fl, fr, fu, fd, theta, U, sumRHO, diffRHO;
double lenx, leny, len, lenxxy, lenxyy;
// Some minor input/output error checking
if (!(nrhs==16)) {
mexErrMsgTxt("Need 16 input arguments");
}
if (!(nlhs==1)) {
mexErrMsgTxt("Need 1 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, dt, 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] );
rho = (double*) mxGetPr( prhs[5] );
u = (double*) mxGetPr( prhs[6] );
v = (double*) mxGetPr( prhs[7] );
NodeP = (int*) mxGetPr( prhs[8] );
Nodeu = (int*) mxGetPr( prhs[9] );
Nodev = (int*) mxGetPr( prhs[10] );
Nbcs = (int) mxGetScalar( prhs[11] );
NbcsDP = (int) mxGetScalar( prhs[12] );
NbcsDu = (int) mxGetScalar( prhs[13] );
NbcsDv = (int) mxGetScalar( prhs[14] );
NBlen = (double) mxGetScalar( prhs[15] );
// 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 rho and u at
// UU
// U
// LL L X R RR
// D
// DD
// with respect to the current cell
#define RHOM rho[NodeP[j + i * (Ny + 2)] - 1]
#define RHOUU rho[NodeP[j - 2 + i * (Ny + 2)] - 1]
#define RHOU rho[NodeP[j - 1 + i *(Ny + 2)] - 1]
#define RHOD rho[NodeP[j + 1 + i * (Ny + 2)] - 1]
#define RHODD rho[NodeP[j + 2 + i * (Ny + 2)] - 1]
#define RHOLL rho[NodeP[j + (i - 2) * (Ny + 2)] - 1]
#define RHOL rho[NodeP[j + (i - 1) * (Ny + 2)] - 1]
#define RHOR rho[NodeP[j + (i + 1) * (Ny + 2)] - 1]
#define RHORR rho[NodeP[j + (i + 2) * (Ny + 2)] - 1]
#define RHORD rho[NodeP[j + 1 + (i + 1) * (Ny + 2)] - 1]
#define RHORU rho[NodeP[j - 1 + (i + 1) * (Ny + 2)] - 1]
#define RHOLD rho[NodeP[j + 1 + (i - 1) * (Ny + 2)] - 1]
#define RHOLU rho[NodeP[j - 1 + (i - 1) * (Ny + 2)] - 1]
#define RHORRD rho[NodeP[j + 1 + (i + 2) * (Ny + 2)] - 1]
#define RHORRU rho[NodeP[j - 1 + (i + 2) * (Ny + 2)] - 1]
#define RHOLLD rho[NodeP[j + 1 + (i - 2) * (Ny + 2)] - 1]
#define RHOLLU rho[NodeP[j - 1 + (i - 2) * (Ny + 2)] - 1]
#define RHORDD rho[NodeP[j + 2 + (i + 1) * (Ny + 2)] - 1]
#define RHORUU rho[NodeP[j - 2 + (i + 1) * (Ny + 2)] - 1]
#define RHOLDD rho[NodeP[j + 2 + (i - 1) * (Ny + 2)] - 1]
#define RHOLUU rho[NodeP[j - 2 + (i - 1) * (Ny + 2)] - 1]
#define RHORRDD rho[NodeP[j + 2 + (i + 2) * (Ny + 2)] - 1]
#define RHORRUU rho[NodeP[j - 2 + (i + 2) * (Ny + 2)] - 1]
#define RHOLLDD rho[NodeP[j + 2 + (i - 2) * (Ny + 2)] - 1]
#define RHOLLUU rho[NodeP[j - 2 + (i - 2) * (Ny + 2)] - 1]
#define UU v[Nodev[j - 1 + (i) * (Ny + 1)] - 1]
#define UD v[Nodev[j + (i) * (Ny + 1)] - 1]
#define UL u[Nodeu[j + (i - 1) * (Ny + 2)] - 1]
#define UR u[Nodeu[j + (i)* (Ny + 2)] - 1]
// ALLOCATE AND POINT TO OUTPUTS
plhs[0] = mxCreateDoubleMatrix((Nx)* Ny , 1, mxREAL);
Fu = (double*)mxGetPr(plhs[0]);
// Calculate lengths to determine narrow band region
lenx = dx*NBlen;
leny = dy*NBlen;
len = sqrt(lenx * lenx + leny * leny);
lenxxy = sqrt(4 * lenx * lenx + leny * leny);
lenxyy = sqrt(lenx * lenx + 4 * leny * leny);
//for(i=1;i<Nx+1;i++) {
// for(j=1;j<Ny+1;j++) {
//TO avoid segfaults, only go within 2 CV's of the boundary
for(i=2; i<Nx; i++) {
for(j=2; j<Ny; j++) {
//If this is a boundary node, don't do any work
if (NodeP[j +(i)*(Ny+2)]<=Nbcs) {
Fu[j-1+(i-1)*Ny] = 0.0;
}
else {
//If this is in the narrowband do work
double lenx;
double leny;
lenx = dx * NBlen;
leny = dy * NBlen;
if ((fabs(RHOM) < MAX(lenx, leny)) || \
(fabs(RHOU) < leny) || (fabs(RHOUU) < 2 * leny) || \
(fabs(RHOD) < leny) || (fabs(RHODD) < 2 * leny) || \
(fabs(RHOL) < lenx) || (fabs(RHOLL) < 2 * lenx) || \
(fabs(RHOR) < lenx) || (fabs(RHORR) < 2 * lenx) || \
(fabs(RHORD) < len) || (fabs(RHORU) < len) || \
(fabs(RHOLD) < len) || (fabs(RHOLU) < len) || \
(fabs(RHORRD) < lenxxy)||(fabs(RHORRU) < lenxxy)|| \
(fabs(RHOLLD) < lenxxy)||(fabs(RHOLLU) < lenxxy)|| \
(fabs(RHORDD) < lenxyy)||(fabs(RHORUU) < lenxyy)|| \
(fabs(RHOLDD) < lenxyy)||(fabs(RHOLUU) < lenxyy)|| \
(fabs(RHORRDD) < 2 * len)||(fabs(RHORRUU) < 2 * len)|| \
(fabs(RHOLLDD) < 2 * len)||(fabs(RHOLLUU) < 2 * len)) {
//do work
//U advection, TVD scheme
/////////////////// WEST or LEFT FACE //////////////////
//mexPrintf("(%d,%d) west\n",i,j);
if(NodeP[j +(i-1)*(Ny+2)]>Nbcs) {
//West face
U = UL;
sumRHO = RHOM + RHOL;
diffRHO = RHOM - RHOL;
theta = ((U + fabs(U)) / 2.0 * (RHOL - RHOLL) \
+ (U - fabs(U)) / 2.0 * (RHOR - RHOM)) \
/ (diffRHO) / U;}
else //treat boundary condition
{
if (Nodeu[j +(i-1)*(Ny+2)]<=NbcsDu || Nodeu[j +(i-1)*(Ny+2)]>Nbcs) {
//West face
U = UL;
}
else {
U = UL + UR;
}
if (NodeP[j +(i-1)*(Ny+2)]<=NbcsDP) {
//West face
sumRHO = RHOM + RHOL;
diffRHO = RHOM - RHOL;
theta = ((U + fabs(U)) / 2.0 * (RHOM - RHOL) \
+ (U - fabs(U)) / 2.0 * (RHOR - RHOM)) \
/ (diffRHO) / U;
}
else { //Neumann BC
sumRHO = 2.0 * RHOM + 0.5 * RHOL;
diffRHO = RHOM - (RHOM + 0.5 * RHOL);
theta = ((U + fabs(U)) / 2.0 * (RHOM - (RHOM + 0.5 * RHOL)) \
+ (U - fabs(U)) / 2.0 * (RHOR - RHOM)) \
/ (diffRHO) / U;
}
}
//Treat additional boundary problems. The theta in TVD
// won't be correct if neumann boundaries are at R
if (i-2>=0) { //This is to avoid segfaults
if(NodeP[j+(i-2)*(Ny+2)]<=Nbcs && NodeP[j+(i-2)*(Ny+2)]>NbcsDP) {
//if the LL boundary is Neumann, use
theta = ((U + fabs(U)) / 2.0 * (RHOL - (RHOL + 0.5 * RHOLL)) \
+ (U - fabs(U)) / 2.0 * (RHOR - RHOM)) \
/ (diffRHO) / U;
}
}
if(NodeP[j+(i+1)*(Ny+2)]<=Nbcs && NodeP[j+(i+1)*(Ny+2)]>NbcsDP) {
//if the most RR boundary is Neumann use:
theta = ((U + fabs(U)) / 2.0 * (RHOL - RHOLL) \
+ (U - fabs(U)) / 2.0 * ((RHOM + 0.5 * RHOR) - RHOM)) \
/ (diffRHO) / U;
}
//TapModification
if(NodeP[j+(i-1)*(Ny+2)]<=Nbcs) //if left cell falls inside boundaries
{theta=0;sumRHO=0;diffRHO=0;}
fl = 0.5 * U * (sumRHO) \
+ 0.5 * fabs(U) * (-1.0 + (1.0 - fabs(U * dt / dx)) * TVD(theta)) \
* (diffRHO);
/////////////////// EAST or RIGHT FACE //////////////////
//mexPrintf("(%d,%d) east\n",i,j);
if(NodeP[j +(i+1)*(Ny+2)]>Nbcs) {
//East face
U = UR;
sumRHO = RHOR + RHOM;
diffRHO = RHOR - RHOM;
theta = ((U + fabs(U)) / 2.0 * (RHOM - RHOL) \
+ (U - fabs(U)) / 2.0 * (RHORR - RHOR)) \
/ (diffRHO) / U;
}
else //treat boundary condition
{
if (Nodeu[j +(i)*(Ny+2)]<=NbcsDu || Nodeu[j +(i)*(Ny+2)]>Nbcs) {
//West face
U = UR;
}
else {
U = UL + UR;
}
if (NodeP[j +(i+1)*(Ny+2)]<=NbcsDP) {
//East face
sumRHO = RHOR + RHOM;
diffRHO = RHOR - RHOM;
theta = ((U + fabs(U)) / 2.0 * (RHOM - RHOL) \
+ (U - fabs(U)) / 2.0 * (RHOR - RHOM)) \
/ (diffRHO) / U;
}
else { //Neumann BC
// mexPrintf("(%d,%d) in sNBC\n",i,j);
sumRHO = 2.0 * RHOM + 0.5 * RHOR;
diffRHO = (RHOM + 0.5 * RHOR) - RHOM;
theta = ((U + fabs(U)) / 2.0 * (RHOM - RHOL) \
+ (U - fabs(U)) / 2.0 * ((RHOM + 0.5 * RHOR) - RHOM)) \
/ (diffRHO) / U;
}
}
//Treat additional boundary problems. The theta in TVD
// won't be correct if neumann boundaries are at
if (i + 2 <= Nx + 1) { //This is to avoid segfaults
if(NodeP[j+(i+2)*(Ny+2)]<=Nbcs && NodeP[j+(i+2)*(Ny+2)]>NbcsDP) {
//if the RR boundary is Neumann, use
theta = ((U + fabs(U)) / 2.0 * (RHOM - RHOL) \
+ (U - fabs(U)) / 2.0 * ((RHOR + 0.5*RHORR) - RHOR)) \
/ (diffRHO) / U;
}
}
if(NodeP[j+(i-1)*(Ny+2)]<=Nbcs && NodeP[j+(i-1)*(Ny+2)]>NbcsDP) {
//if the most L boundary is Neumann use:
theta = ((U + fabs(U)) / 2.0 * (RHOM - (RHOM + 0.5 * RHOL)) \
+ (U - fabs(U)) / 2.0 * (RHORR - RHOR)) \
/ (diffRHO) / U;
}
//TapModification
if(NodeP[j+(i+1)*(Ny+2)]<=Nbcs) //if right cell falls inside boundaries
{theta=0;sumRHO=0;diffRHO=0;}
fr = 0.5 * U * (sumRHO) \
+ 0.5 * fabs(U) * (-1.0 + (1.0 - fabs(U * dt / dx)) * TVD(theta)) \
* (diffRHO);
/////////////////// SOUTH or DOWN FACE //////////////////
if (NodeP[j+1+ i*(Ny+2)]>Nbcs) {
//South Face
U = UD;
sumRHO = RHOD + RHOM;
diffRHO = RHOM - RHOD;
theta = ((U + fabs(U)) / 2.0 * (RHOD - RHODD) \
+ (U - fabs(U)) / 2.0 * (RHOU - RHOM)) \
/ (diffRHO) / U;
}
else { //treat boundaries
if (Nodev[j +(i)*(Ny+1)]<=NbcsDv || Nodev[j +(i)*(Ny+1)]>Nbcs) { //Dirichlet on v
U = UD;
}
else {
U = UD + UU;
}
if (NodeP[j+1+(i)*(Ny+2)]<=NbcsDP) { //Dirichlet on v
sumRHO = RHOD + RHOM;
diffRHO = RHOM - RHOD;
theta = ((U + fabs(U)) / 2.0* (RHOM - RHOD) \
+ (U - fabs(U)) / 2.0 * (RHOU - RHOM)) \
/ (diffRHO) / U;
}
else { //Neumann on u
sumRHO = RHOD + 2.0 * RHOM;
diffRHO = RHOM - (RHOM + 0.5 * RHOD);
theta = ((U + fabs(U)) / 2.0 * (RHOM - (RHOM + 0.5 * RHOD)) \
+ (U - fabs(U)) / 2.0 * (RHOU - RHOM)) \
/ (diffRHO) / U;
}
}
//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(NodeP[j+2+i*(Ny+2)]<=Nbcs && NodeP[j+2+i*(Ny+2)]>NbcsDP) {
//if the most south boundary is Neumann use:
theta = ((U + fabs(U)) / 2.0 * (RHOD - (RHOD + 0.5 * RHODD)) \
+ (U - fabs(U)) / 2.0 * (RHOU - RHOM)) \
/ (diffRHO) / U;
}
}
if(NodeP[j-1+i*(Ny+2)]<=Nbcs && NodeP[j-1+i*(Ny+2)]>NbcsDP) {
//if the most north boundary is Neumann use:
theta = ((U + fabs(U)) / 2.0 * (RHOD - RHODD) \
+ (U - fabs(U)) / 2.0 * ((RHOM + 0.5 * RHOM) - RHOM)) \
/ (diffRHO) / U;
}
//TapModification
if(NodeP[j+1+ i*(Ny+2)]<=Nbcs) //if south cell falls inside boundaries
{theta=0;sumRHO=0;diffRHO=0;}
fd = 0.5 * U*(sumRHO) \
+ 0.5 * fabs(U) * (-1.0 + (1.0 - fabs(U * dt / dy)) * TVD(theta)) \
* (diffRHO);
/////////////////// NORTH or UP FACE //////////////////
//North Face
if (NodeP[j-1+ i*(Ny+2)]>Nbcs) {
U = UU;
sumRHO = RHOU + RHOM;
diffRHO = RHOU - RHOM;
theta = ((U + fabs(U)) / 2.0 * (RHOM - RHOD) \
+ (U - fabs(U)) / 2.0 * (RHOUU - RHOU)) \
/ (diffRHO) / U;
}
else { //treat boundaries
if (Nodev[j-1 +(i)*(Ny+1)]<=NbcsDv || Nodev[j-1 +(i)*(Ny+1)]>Nbcs) { //Dirichlet on v
U = UU;
}
else { //Neuman on v
U = UU + UD;
}
if (NodeP[j-1+(i)*(Ny+2)]<=NbcsDP) { //Dirichlet on u
sumRHO = RHOU + RHOM;
diffRHO = RHOU - RHOM;
theta = ((U + fabs(U)) / 2.0 * (RHOM - RHOD) \
+ (U - fabs(U)) / 2.0 * (RHOU - RHOM)) \
/ (diffRHO) / U;
}
else { //Neuman on u
sumRHO = 0.5 * RHOU + 2.0 * RHOM;
diffRHO = (0.5 * RHOU + RHOM) - RHOM;
theta = ((U + fabs(U)) / 2.0 * (RHOM - RHOD) \
+ (U - fabs(U)) / 2.0 * ((0.5 * RHOU + RHOM) - RHOM)) \
/ (diffRHO) / U;
}
}
//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(NodeP[j-2+i*(Ny+2)]<=Nbcs && NodeP[j-2+i*(Ny+2)]>NbcsDP) {
//if the most north boundary is Neumann use
theta = ((U + fabs(U)) / 2.0 * (RHOM - RHOD) \
+ (U - fabs(U)) / 2.0 * ((RHOU + 0.5 * RHOUU) - RHOU)) \
/ (diffRHO) / U;
}
}
if(NodeP[j+1+i*(Ny+2)]<=Nbcs && NodeP[j+1+i*(Ny+2)]>NbcsDP) {
//if the most south boundary is Neumann use
theta = ((U + fabs(U)) / 2.0 * (RHOM - (RHOM + 0.5 * RHOD)) \
+ (U - fabs(U)) / 2.0 * (RHOUU - RHOU)) \
/ (diffRHO) / U;
}
//TapModification
if(NodeP[j-1+ i*(Ny+2)]<=Nbcs) //if left cell falls inside boundaries
{theta=0;sumRHO=0;diffRHO=0;}
fu = 0.5 * U*(sumRHO) \
+ 0.5 * fabs(U) * (-1.0 + (1.0 - fabs(U * dt / dy)) * TVD(theta)) \
* (diffRHO);
//Calculate Total flux change for this cell
Fu[j-1+(i-1)*Ny] = (fl - fr) / dx + (fd - fu) / dy;
}
else {
Fu[j-1+(i-1)*Ny] = 0.0;
}
}
}
}
}