start lbm version

Makefile

+.PHONY: all clean
+
+DEBUG    := n
+SANITIZE := n
+
+ifeq ($(SANITIZE), y)
+SANITIZEFLAGS := -fsanitize=address -fsanitize=undefined
+endif
+
+ifeq ($(DEBUG), y)
+WARNFLAGS := -Wall -Wpedantic -Wextra -Waddress -Waggressive-loop-optimizations \
+  -Wcast-qual -Wcast-align -Wmissing-declarations \
+  -Wdouble-promotion -Wuninitialized -Winit-self \
+  -Wstrict-aliasing -Wsuggest-attribute=const -Wtrampolines -Wfloat-equal \
+  -Wshadow -Wunsafe-loop-optimizations -Wfloat-conversion -Wlogical-op \
+  -Wnormalized -Wdisabled-optimization -Whsa -Wconversion -Wunused-result
+
+CPPFLAGS += -DDEBUG
+OPTIFLAGS := -Og
+else
+CPPFLAGS += -DNDEBUG
+OPTIFLAGS := -O3 -march=native -flto -fno-math-errno
+endif
+
+CPPFLAGS += -I include
+CFLAGS += -std=c11 -g $(WARNFLAGS) $(OPTIFLAGS) $(SANITIZEFLAGS) -fopenmp
+LDFLAGS += -g $(OPTIFLAGS) $(SANITIZEFLAGS) -lm -fopenmp
+
+all: transport_openmp
+
+transport_openmp: transport_openmp.o
+	$(CC) $^ $(LDFLAGS) -o $@
+
+transport_openmp.o: transport_openmp.c
+
+clean:
+	rm transport_openmp.o transport_openmp

lbm_cl.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+
+# resolution of a transport equation by the finite volume method
+# on regular grid
+
+# regular python implementation compared to a pyopencl version
+
+
+from __future__ import absolute_import, print_function
+import pyopencl as cl
+import numpy as np
+import matplotlib.pyplot as plt
+
+import time
+
+
+##################" definition of default values
+# number of conservative variables
+_m = 1
+
+# number of kinetic variables
+_n = 4 * _m
+
+
+# grid size
+_nx = 64
+_ny = 64
+
+Lx = 1.
+Ly = 1.
+
+_dx = Lx / _nx
+_dy = Ly / _ny
+
+# transport velocity
+vel = np.array([1., 1.])
+
+# lattice speed
+_vmax = 2.
+
+
+# time stepping
+_Tmax = 0.5 / _vmax
+cfl = 1
+
+_dt = cfl * _dx / _vmax
+############# end of default values
+
+def exact_sol(xy, t):
+    x = xy[0] - t * vel[0] - 0.5
+    y = xy[1] - t * vel[1] - 0.5
+    d2 = x * x + y *y
+    w = np.exp(-30*d2)
+    return w
+
+
+
+def solve_ocl(m = _m, n = _n, nx = _nx, ny = _ny,
+              Tmax = _Tmax, vmax = _vmax,
+              dx = _dx, dy = _dy,
+              dt = _dt, exact_sol = exact_sol,
+              animate = False):
+
+    # load and adjust  C program
+    source = open("lbm_kernels.cl", "r").read()
+    source = source.replace("_nx_", "("+str(nx)+")")
+    source = source.replace("_ny_", "("+str(ny)+")")
+    source = source.replace("_dx_", "("+str(dx)+"f)")
+    source = source.replace("_dy_", "("+str(dy)+"f)")
+    source = source.replace("_dt_", "("+str(dt)+"f)")
+    source = source.replace("_m_", "("+str(m)+")")
+    source = source.replace("_n_", "("+str(n)+")")    
+    source = source.replace("_vx_", "("+str(vel[0])+"f)")
+    source = source.replace("_vy_", "("+str(vel[1])+"f)")
+    source = source.replace("_lambda_", "("+str(vmax)+"f)")
+
+    #print(source)
+
+    #exit(0)
+
+    # OpenCL init
+    ctx = cl.create_some_context()
+    mf = cl.mem_flags
+
+    # compile OpenCL C program
+    prg = cl.Program(ctx, source).build(options = "-cl-strict-aliasing  \
+                                                  -cl-fast-relaxed-math")
+
+    # create OpenCL buffers
+    fn_gpu   = cl.Buffer(ctx, mf.READ_WRITE,
+                         size=(4 * m * nx * ny * np.dtype('float32').itemsize))
+    fnp1_gpu = cl.Buffer(ctx, mf.READ_WRITE,
+                         size=(4 * m * nx * ny * np.dtype('float32').itemsize))
+
+    # create a queue (for submitting opencl operations)
+    queue = cl.CommandQueue(ctx, properties=cl.command_queue_properties.PROFILING_ENABLE)
+
+
+    # init data
+    event = prg.init_sol(queue, (nx * ny, ), (32, ), fn_gpu)
+    event.wait()
+
+    # number of animation frames
+    nbplots = 10
+    itermax = int(np.floor(Tmax / dt))
+    iterplot = int(itermax / nbplots)
+
+    # time loop
+    t = 0
+    iter = 0
+    elapsed = 0.;
+    fn_cpu = np.empty((4 * m * nx * ny, ), dtype = np.float32)
+
+    print("start OpenCL computations...")
+    while t < Tmax:
+        t = t + dt
+        iter = iter + 1
+        #event = prg.time_step(queue, (nx * ny, ), (32, ), wn_gpu, wnp1_gpu)
+        event = prg.time_step(queue, (nx * ny, ), (64, ), fn_gpu, fnp1_gpu)
+        #event = prg.time_step(queue, (nx * ny, ), (32, ), wn_gpu, wnp1_gpu, wait_for = [event])
+        event.wait()
+        elapsed += 1e-9 * (event.profile.end - event.profile.start)
+        # exchange buffer references for avoiding a copy
+        fn_gpu, fnp1_gpu = fnp1_gpu, fn_gpu
+        print("iter=",iter, " t=",t, "elapsed (s)=",elapsed)
+        if iter % iterplot == 0 and animate:
+            cl.enqueue_copy(queue, fn_cpu, fn_gpu).wait()
+            wplot = np.reshape(fn_cpu, (nx, ny))
+            plt.clf()
+            plt.imshow(wplot,vmin=0, vmax=1)
+            plt.gca().invert_yaxis()
+            plt.colorbar()
+            plt.pause(0.01)
+
+    # copy OpenCL data to CPU and return the results
+    cl.enqueue_copy(queue, fn_cpu, fn_gpu).wait()
+
+    wplot_gpu = np.reshape(fn_cpu,(nx, ny))
+    return wplot_gpu
+
+# gpu solve
+wplot_gpu = solve_ocl()
+plt.clf()
+plt.imshow(wplot_gpu, vmin=0, vmax=1)
+plt.gca().invert_yaxis()
+plt.colorbar()
+plt.show()
+
+
+# check difference
+# plt.clf()
+# plt.imshow(wplot_cpu-wplot_gpu)
+# plt.gca().invert_yaxis()
+# plt.colorbar()
+# plt.show()
+
+

lbm_kernels.cl

+#define _NX _nx_
+#define _NY _ny_
+#define _DX _dx_
+#define _DY _dy_
+#define _DT _dt_
+#define _M _m_
+#define _N _n_
+#define _VX _vx_
+#define _VY _vy_
+
+#define _LAMBDA _lambda_
+
+
+#define _VOL (_DX * _DY)
+
+__constant int dir[4][2] = { {-1, 0}, {1, 0},
+			     {0, -1}, {0, 1}};
+
+__constant float ds[4] = { _DY, _DY, _DX, _DX };
+
+
+
+// physical flux of the hyperbolic system
+void flux_phy(const float *w, const float* vnorm, float* flux){
+  float vdotn = _VX * vnorm[0] + _VY * vnorm[1];
+  for(int iv = 0; iv < _M; iv++){
+    flux[iv] = vdotn * w[iv];
+  }
+}
+
+// equilibrium "maxwellian" from macro data w
+void w2f(const float *w, float *f){
+  for(int d = 0; d < 4; d++){
+    float flux[_M];
+    float vnorm[2] = {dir[d][0], dir[d][1]};
+    flux_phy(w, vnorm, flux);   
+    for(int iv = 0; iv < _M; iv++){
+      f[d * _M + iv] = w[iv] + flux[iv] / _LAMBDA;
+    }
+  }
+}
+
+
+// macro data w from micro data f
+void f2w(const float *f, float *w){
+  for(int iv = 0; iv < _M; iv++) w[iv] = 0;
+  for(int d = 0; d < 4; d++){
+    for(int iv = 0; iv < _M; iv++){
+      w[iv] = f[d * _M + iv];
+    }
+  }
+}
+
+
+// exact macro data w
+void exact_sol(float* xy, float t, float* w){
+  float x = xy[0] - t * _VX - 0.5f;
+  float y = xy[1] - t * _VY - 0.5f;
+  float d2 = x * x + y *y;
+  for(int iv = 0; iv < _M; iv++)
+    w[iv] = exp(-30*d2);
+}
+
+// initial condition on the macro data
+__kernel void init_sol(__global  float *fn){
+
+  int id = get_global_id(0);
+  
+  int i = id % _NX;
+  int j = id / _NX;
+
+  int ngrid = _NX * _NY;
+
+  float wnow[_M];
+
+  float t = 0;
+  float xy[2] = {i * _DX + _DX / 2, j * _DY + _DY / 2};
+  
+  exact_sol(xy, t, wnow);
+
+  float fnow[_N];
+  w2f(wnow, fnow);
+
+  //printf("x=%f, y=%f \n",xy[0],xy[1]);
+  // load middle value
+  for(int ik = 0; ik < _N; ik++){
+    int imem = i + j * _NX + ik * ngrid;
+    fn[imem] =  fnow[ik];
+  }
+
+}
+
+
+
+// one time step of the LBM scheme
+__kernel void time_step(__global  float *fn, __global float *fnp1){
+
+  int id = get_global_id(0);
+  
+  int i = id % _NX;
+  int j = id / _NX;
+
+
+  int ngrid = _NX * _NY;
+
+  float fnow[_N];
+  
+  // periodic shift in each direction
+  for(int d = 0; d < 4; d++){
+    int iR = (i - dir[d][0]) % _NX;
+    int jR = (j - dir[d][1]) % _NY;
+    
+    for(int iv = 0; iv < _M; iv++){
+      int ik = d * _M + iv;
+      int imem = iR + jR * _NX + ik * ngrid;
+      fnow[ik] = fn[imem];
+    }
+  }
+
+  float wnow[_M];
+  f2w(fnow, wnow);
+  
+  float fnext[_N];
+  // first order relaxation
+  w2f(wnow, fnext);
+
+  // second order relaxation
+  for(int ik = 0; ik < _N; ik++)
+    fnext[ik] = 2 * fnext[ik] - fnow[ik];
+
+  // save
+  for(int ik = 0; ik < _M; ik++){
+      int imem = i + j * _NX + ik * ngrid;
+      fnp1[imem] = fnext[ik];
+  }
+
+
+  
+
+}

prog_lukas/Equations.cpp

+/*
+ * Equations.cpp
+ *
+ *  Created on: Oct 4, 2017
+ *      Author: lukas
+ */
+
+#include "Equations.h"
+
+Equations::Equations(int NumConVar,double Lambda,double* VN1,double* VN2) {
+
+	numberOfConsVar = NumConVar;
+	normalVector1 = VN1;
+	normalVector2 = VN2;
+	lambda = Lambda;
+
+}
+
+double* Equations::Cons2EquilibriumFunc(Grid* grid,int cellIndexI,int cellIndexJ){
+
+	double* equilibriumFunctions = new double[numberOfConsVar * grid->numberOfdirections];
+	double* flux1 = FluxInDirection1(grid,cellIndexI,cellIndexJ);
+
+    if(grid->numberOfdirections == 5){
+    	double* flux2 = FluxInDirection2(grid,cellIndexI,cellIndexJ);
+		for(int c = 0; c < numberOfConsVar; c++){
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionRight] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 5.0 - flux1[c] / 2.0 / lambda;
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionLeft] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 5.0 + flux1[c] / 2.0 / lambda;
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionUp] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 5.0 - flux2[c] / 2.0 / lambda;
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionDown] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 5.0 + flux2[c] / 2.0 / lambda;
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionDown + 1] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 5.0;
+		}
+		delete[] flux2;
+		flux2 = 0;
+    }else if(grid->numberOfdirections == 4){
+    	double* flux2 = FluxInDirection2(grid,cellIndexI,cellIndexJ);
+		for(int c = 0; c < numberOfConsVar; c++){
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionRight] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 4.0 - flux1[c] / 2.0 / lambda;
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionLeft] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 4.0 + flux1[c] / 2.0 / lambda;
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionUp] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 4.0 - flux2[c] / 2.0 / lambda;
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionDown] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 4.0 + flux2[c] / 2.0 / lambda;
+		}
+		delete[] flux2;
+		flux2 = 0;
+    }else if(grid->numberOfdirections == 3){
+    	for(int c = 0; c < numberOfConsVar; c++){
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionRight] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 3.0 - flux1[c] / 2.0 / lambda;
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionLeft] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 3.0 + flux1[c] / 2.0 / lambda;
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionLeft + 1] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 3.0;
+    	}
+    }
+    else if(grid->numberOfdirections == 2){
+		for(int c = 0; c < numberOfConsVar; c++){
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionRight] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 2.0 - flux1[c] / 2.0 / lambda;
+			equilibriumFunctions[c * grid->numberOfdirections + grid->indexDirectionLeft] = grid->conservativeVariables[cellIndexI][cellIndexJ][c] / 2.0 + flux1[c] / 2.0 / lambda;
+		}
+	}
+    else{
+    	printf("Number of directions is neither 2, 3 nor 4! -> check directions/dimensions \n");
+    	exit(EXIT_FAILURE);
+    }
+	delete[] flux1;
+	flux1 = 0;
+
+	return equilibriumFunctions;
+}
+
+double* Equations::FluxInDirection1(Grid* grid,int cellIndexI,int cellIndexJ){return Flux(grid,normalVector1,cellIndexI,cellIndexJ);}
+double* Equations::FluxInDirection2(Grid* grid,int cellIndexI,int cellIndexJ){return Flux(grid,normalVector2,cellIndexI,cellIndexJ);}
+
+
+void Equations::computeInitialKineticData(Grid* grid){
+
+	for(int i = 0; i < grid->Nx; i++){
+			for(int j = 0; j < grid->Ny; j++){
+				grid->kineticVariables[i][j] = Cons2EquilibriumFunc(grid,i,j);
+			}
+	}
+
+}

prog_lukas/Equations.h

+/*
+ * Equations.h
+ *
+ *  Created on: Oct 4, 2017
+ *      Author: lukas
+ */
+
+#ifndef EQUATIONS_H_
+#define EQUATIONS_H_
+
+#include "Grid.h"
+
+class Equations {
+
+	public:
+		Equations(int NumConVar,double Lambda,double* VN1,double* VN2);
+
+		int numberOfConsVar;
+		double* normalVector1;
+		double* normalVector2;
+		double lambda;
+
+		virtual double* Flux(Grid* Grid, double* NormalVector, int CellIndexI, int CellIndexJ) =0;
+		double* FluxInDirection1(Grid* Grid, int CellIndexI, int CellIndexJ);
+		double* FluxInDirection2(Grid* Grid, int CellIndexI, int CellIndexJ);
+
+		double* Cons2EquilibriumFunc(Grid* Grid, int CellIndexI, int CellIndexJ);
+		void computeInitialKineticData(Grid* grid);
+		virtual ~Equations() {};
+
+
+};
+
+#endif /* EQUATIONS_H_ */

prog_lukas/EquationsMHD.cpp

+#include "EquationsMHD.h"
+
+EquationsMHD::EquationsMHD(int NumConVar,double Lambda,double DCS,double G,double* VN1,double* VN2):Equations(NumConVar,Lambda,VN1,VN2) {
+
+	divergenceCleaningSpeed =DCS;
+	Gamma = G;
+
+	std::cout<<"Equations successfully created"<<std::endl;
+	std::cout<<"-> Type: Magnetohydrodynamics (MHD)"<<std::endl;
+
+}
+
+double* EquationsMHD::Flux(Grid* grid,double* normalVector,int cellIndexI,int cellIndexJ){
+
+	double* flux = new double[numberOfConsVar];
+
+	double un = grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumX] / grid->conservativeVariables[cellIndexI][cellIndexJ][indexDensity] * normalVector[grid->indexDimensionX] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumY] / grid->conservativeVariables[cellIndexI][cellIndexJ][indexDensity] * normalVector[grid->indexDimensionY] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumZ] / grid->conservativeVariables[cellIndexI][cellIndexJ][indexDensity] * normalVector[grid->indexDimensionZ];
+
+	double bn = grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] * normalVector[grid->indexDimensionX] + grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] * normalVector[grid->indexDimensionY] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ] * normalVector[grid->indexDimensionZ];
+
+	double p = (Gamma - 1.0) *(grid->conservativeVariables[cellIndexI][cellIndexJ][indexEnergy] -  (grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumX] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumX] / grid->conservativeVariables[cellIndexI][cellIndexJ][indexDensity] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumY] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumY] / grid->conservativeVariables[cellIndexI][cellIndexJ][indexDensity] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumZ] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumZ] / grid->conservativeVariables[cellIndexI][cellIndexJ][indexDensity]) / 2.0 -
+			(grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] +
+	   		grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ]) / 2.0);
+
+	flux[indexDensity] = grid->conservativeVariables[cellIndexI][cellIndexJ][indexDensity] * un;
+
+	flux[indexMomentumX] = un * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumX] +
+			(p + (grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ]) / 2.0) * normalVector[grid->indexDimensionX] -
+			bn * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX];
+
+	flux[indexMomentumY] = un * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumY] +
+			(p + (grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ]) / 2.0) * normalVector[grid->indexDimensionY] -
+			bn * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY];
+
+	flux[indexMomentumZ] = un * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumZ] +
+			(p + (grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ]) / 2.0) * normalVector[grid->indexDimensionZ] -
+			bn * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ];
+
+
+	flux[indexEnergy] = (grid->conservativeVariables[cellIndexI][cellIndexJ][indexEnergy] + p + (grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ]) / 2.0) * un -
+			(grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumX] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumY] +
+			grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ] * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumZ]) / grid->conservativeVariables[cellIndexI][cellIndexJ][indexDensity] * bn;
+
+	flux[indexMagneticFieldX] = -bn * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumX] / grid->conservativeVariables[cellIndexI][cellIndexJ][indexDensity] +
+			un * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldX] + grid->conservativeVariables[cellIndexI][cellIndexJ][indexDivergenceCleaningPotential] * normalVector[grid->indexDimensionX];
+
+	flux[indexMagneticFieldY] = -bn * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumY] / grid->conservativeVariables[cellIndexI][cellIndexJ][indexDensity] +
+			un * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldY] + grid->conservativeVariables[cellIndexI][cellIndexJ][indexDivergenceCleaningPotential] * normalVector[grid->indexDimensionY];
+
+	flux[indexMagneticFieldZ] = -bn * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMomentumZ] / grid->conservativeVariables[cellIndexI][cellIndexJ][indexDensity] +
+			un * grid->conservativeVariables[cellIndexI][cellIndexJ][indexMagneticFieldZ] + grid->conservativeVariables[cellIndexI][cellIndexJ][indexDivergenceCleaningPotential] * normalVector[grid->indexDimensionZ];
+
+	flux[indexDivergenceCleaningPotential] = divergenceCleaningSpeed * divergenceCleaningSpeed * bn;
+
+	return flux;
+}
+
+
+double EquationsMHD::computeAbsoluteSpeedSquared(double* conservativeVariables){
+
+	double u = conservativeVariables[indexMomentumX]/conservativeVariables[indexDensity];
+	double v = conservativeVariables[indexMomentumY]/conservativeVariables[indexDensity];
+	double w = conservativeVariables[indexMomentumZ]/conservativeVariables[indexDensity];
+
+	return u*u+v*v+w*w;
+}
+
+double EquationsMHD::computeAbsoluteMagneitcFieldSquared(double* conservativeVariables){
+
+	return	conservativeVariables[indexMagneticFieldX]*conservativeVariables[indexMagneticFieldX]+
+			conservativeVariables[indexMagneticFieldY]*conservativeVariables[indexMagneticFieldY]+
+			conservativeVariables[indexMagneticFieldZ]*conservativeVariables[indexMagneticFieldZ];
+
+}
+
+double EquationsMHD::computeAbsoluteSpeed(double* conservativeVariables){return sqrt(computeAbsoluteSpeedSquared(conservativeVariables));}
+double EquationsMHD::computeAbsoluteMagneticField(double* conservativeVariables){return sqrt(computeAbsoluteMagneitcFieldSquared(conservativeVariables));}
+
+double EquationsMHD::computeAbsoluteMomentum(double* conservativeVariables){return computeAbsoluteSpeed(conservativeVariables) * conservativeVariables[indexDensity];}
+
+double EquationsMHD::computePressure(double* conservativeVariables){
+
+	double u = conservativeVariables[indexMomentumX]/conservativeVariables[indexDensity];
+	double v = conservativeVariables[indexMomentumY]/conservativeVariables[indexDensity];
+	double w = conservativeVariables[indexMomentumZ]/conservativeVariables[indexDensity];
+	double absoluteSpeedSquared = u*u+v*v+w*w;
+	double absoluteMagneticFieldSquared = conservativeVariables[indexMagneticFieldX]*conservativeVariables[indexMagneticFieldX]+conservativeVariables[indexMagneticFieldY]*conservativeVariables[indexMagneticFieldY]+conservativeVariables[indexMagneticFieldZ]*conservativeVariables[indexMagneticFieldZ];
+
+	return (conservativeVariables[indexEnergy] - (conservativeVariables[indexDensity] * absoluteSpeedSquared/2.0)-(absoluteMagneticFieldSquared/2.0))  *(Gamma -1.0);
+}
+
+double EquationsMHD::computeSpeedOfSound(double* conservativeVariables){return sqrt(Gamma * computePressure(conservativeVariables) / conservativeVariables[indexDensity]);}
+double EquationsMHD::computeMach(double* conservativeVariables){return computeAbsoluteSpeed(conservativeVariables)/computeSpeedOfSound(conservativeVariables);}
+
+double EquationsMHD::computeCharacteristicVariableEuler1(double* conservativeVariables){return computePressure(conservativeVariables)+ computeAbsoluteMomentum(conservativeVariables) * computeSpeedOfSound(conservativeVariables);}
+double EquationsMHD::computeCharacteristicVariableEuler2(double* conservativeVariables){return computePressure(conservativeVariables)- computeAbsoluteMomentum(conservativeVariables) * computeSpeedOfSound(conservativeVariables);}
+double EquationsMHD::computeCharacteristicVariableEuler3(double* conservativeVariables){return computePressure(conservativeVariables)- conservativeVariables[indexDensity] * computeSpeedOfSound(conservativeVariables) * computeSpeedOfSound(conservativeVariables);}
+
+double* EquationsMHD::computeCharacteristicVariablesEuler(double* conservativeVariables){
+
+	double* characteristicVariables = new double[3];
+
+	characteristicVariables[0] = computeCharacteristicVariableEuler1(conservativeVariables);
+	characteristicVariables[1] = computeCharacteristicVariableEuler2(conservativeVariables);
+	characteristicVariables[2] = computeCharacteristicVariableEuler3(conservativeVariables);
+
+	return characteristicVariables;
+}
+
+double EquationsMHD::computeEulerEVi(double* conservativeVariables,int MomentumIndex,int i){
+	double eigenvalue = 0.0;
+	if(i == 1){eigenvalue = computeEulerEV1(conservativeVariables,MomentumIndex);}
+	if(i == 2){eigenvalue = computeEulerEV2(conservativeVariables,MomentumIndex);}
+	if(i == 3){eigenvalue = computeEulerEV3(conservativeVariables,MomentumIndex);}
+	return eigenvalue;
+}
+
+double EquationsMHD::computeEulerEV1(double* conservativeVariables,int MomentumIndex){
+