From aa66aadb7451fe3f82f4c05e2fe860084e1e644f Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian@jauslin.org>
Date: Thu, 26 May 2022 15:16:44 -0400
Subject: Driving force as vector instead of function

---
 src/driving.c       | 50 +++++++++++++++++++++++++++++++-------------------
 src/driving.h       |  6 ++----
 src/main.c          | 41 +++++++++++++++++++++++++++--------------
 src/navier-stokes.c | 18 +++++++++---------
 src/navier-stokes.h | 14 +++++++-------
 5 files changed, 76 insertions(+), 53 deletions(-)

(limited to 'src')

diff --git a/src/driving.c b/src/driving.c
index f12c5fe..70c2564 100644
--- a/src/driving.c
+++ b/src/driving.c
@@ -1,28 +1,40 @@
 #include "driving.h"
+#include "navier-stokes.h"
 #include <math.h>
 
-_Complex double g_test(
-  int kx,
-  int ky
+// test driving function
+int g_test(
+  _Complex double* g,
+  int K1,
+  int K2
 ){
-  //return sqrt(kx*kx*ky*ky)*exp(-(kx*kx+ky*ky));
-  if(kx==2 && ky==-1){
-    return 0.5+sqrt(3)/2*I;
+  int kx,ky;
+  for(kx=-K1;kx<=K1;kx++){
+    for (ky=-K2;ky<=K2;ky++){
+      if(kx==2 && ky==-1){
+	g[klookup(kx,ky,K1,K2)]=0.5+sqrt(3)/2*I;
+      }
+      else if(kx==-2 && ky==1){
+        g[klookup(kx,ky,K1,K2)]=0.5-sqrt(3)/2*I;
+      }
+      else{
+        g[klookup(kx,ky,K1,K2)]=0.;
+      }
+    }
   }
-  else if(kx==-2 && ky==1){
-     return 0.5-sqrt(3)/2*I;
-  }
-  return 0.;
+
+  return 0;
 }
 
-#define UNUSED(x) (void)(x)
-_Complex double g_zero(
-  int kx,
-  int ky
+int g_zero(
+  _Complex double* g,
+  int K1,
+  int K2
 ){
-  // avoid unused variable warnings
-  UNUSED(kx);
-  UNUSED(ky);
-  return 0.;
-}
+  int i;
+  for(i=0;i<(2*K1+1)*(2*K2+1);i++){
+    g[i]=0.;
+  }
 
+  return 0;
+}
diff --git a/src/driving.h b/src/driving.h
index 5d3b415..1d07812 100644
--- a/src/driving.h
+++ b/src/driving.h
@@ -1,9 +1,7 @@
 #ifndef DRIVING_H
 #define DRIVING_H
 
-#include <complex.h>
-
-_Complex double g_zero( int kx, int ky);
-_Complex double g_test( int kx, int ky);
+int g_zero(_Complex double* g, int K1, int K2);
+int g_test(_Complex double* g, int K1, int K2);
 
 #endif
diff --git a/src/main.c b/src/main.c
index 27ed314..f074720 100644
--- a/src/main.c
+++ b/src/main.c
@@ -32,6 +32,8 @@ int read_args(int argc, const char* argv[], char** params, unsigned int* driving
 int read_params(char* param_str, nstrophy_parameters* parameters);
 int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2);
 
+// set driving force
+_Complex double* set_driving(unsigned int driving, nstrophy_parameters parameters);
 // set initial condition
 _Complex double* set_init(unsigned int init, nstrophy_parameters parameters);
 
@@ -53,11 +55,11 @@ int main (
 ){
   char* param_str=NULL;
   nstrophy_parameters parameters;
-  _Complex double (*g)(int,int);
   int ret;
   unsigned int driving,command,init;
   unsigned int nthreads=1;
   _Complex double* u0;
+  _Complex double *g;
 
   command=0;
   driving=0;
@@ -76,19 +78,7 @@ int main (
   }
 
   // set driving force
-  switch(driving){
-  case DRIVING_ZERO:
-    g=g_zero;
-    break;
-  case DRIVING_TEST:
-    g=g_test;
-    break;
-
-  default:
-    g=g_zero;
-    break;
-  }
-
+  g=set_driving(driving, parameters);
   // set initial condition
   u0=set_init(init, parameters);
 
@@ -110,6 +100,7 @@ int main (
     print_usage();
   }
 
+  free(g);
   free(u0);
 
   return(0);
@@ -441,6 +432,28 @@ int set_parameter(
   return(0);
 }
 
+// set driving force
+_Complex double* set_driving(
+  unsigned int driving,
+  nstrophy_parameters parameters
+){
+  _Complex double* g=calloc(sizeof(_Complex double),(2*parameters.K1+1)*(2*parameters.K2+1));
+
+  switch(driving){
+  case DRIVING_ZERO:
+    g_zero(g, parameters.K1, parameters.K2);
+    break;
+  case DRIVING_TEST:
+    g_test(g, parameters.K1, parameters.K2);
+    break;
+
+  default:
+    g_test(g, parameters.K1, parameters.K2);
+    break;
+  }
+   return g;
+}
+
 // set initial condition
 _Complex double* set_init(
   unsigned int init,
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index c52b810..fadc87e 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -13,7 +13,7 @@ int uk(
   double delta,
   double L,
   _Complex double* u0,
-  _Complex double (*g)(int,int),
+  _Complex double* g,
   unsigned int print_freq,
   unsigned int nthreads
 ){
@@ -80,7 +80,7 @@ int energy(
   double delta,
   double L,
   _Complex double* u0,
-  _Complex double (*g)(int,int),
+  _Complex double* g,
   unsigned int print_freq,
   unsigned int nthreads
 ){
@@ -131,7 +131,7 @@ int enstrophy(
   double delta,
   double L,
   _Complex double* u0,
-  _Complex double (*g)(int,int),
+  _Complex double* g,
   unsigned int print_freq,
   unsigned int nthreads
 ){
@@ -185,7 +185,7 @@ int quiet(
   double delta,
   double L,
   _Complex double* u0,
-  _Complex double (*g)(int,int),
+  _Complex double* g,
   unsigned int nthreads
 ){
   _Complex double* u;
@@ -306,7 +306,7 @@ int ins_step(
   double nu,
   double delta,
   double L,
-  _Complex double (*g)(int,int),
+  _Complex double* g,
   fft_vect fft1,
   fft_vect fft2,
   fft_vect ifft,
@@ -383,7 +383,7 @@ int ins_rhs(
   int N2,
   double nu,
   double L,
-  _Complex double (*g)(int,int),
+  _Complex double* g,
   fft_vect fft1,
   fft_vect fft2,
   fft_vect ifft
@@ -413,7 +413,7 @@ int ins_rhs(
     for(ky=-K2;ky<=K2;ky++){
       if(kx!=0 || ky!=0){
 	// enforce the reality of u by adding ifft.fft(k) and the conjugate of ifft.fft(-k)
-        out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*(ifft.fft[klookup(kx,ky,N1,N2)]+conj(ifft.fft[klookup(-kx,-ky,N1,N2)]))/2;
+        out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+g[klookup(kx,ky,2*K1+1,2*K2+1)]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*(ifft.fft[klookup(kx,ky,N1,N2)]+conj(ifft.fft[klookup(-kx,-ky,N1,N2)]))/2;
       }
     }
   }
@@ -535,7 +535,7 @@ _Complex double compute_alpha(
   _Complex double* u,
   int K1,
   int K2,
-  _Complex double (*g)(int,int)
+  _Complex double* g
 ){
   _Complex double num=0;
   _Complex double denom=0;
@@ -544,7 +544,7 @@ _Complex double compute_alpha(
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
       denom+=(kx*kx+ky*ky)*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)])*(1+(ky!=0?kx*kx/ky/ky:0));
-      num+=(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj((*g)(kx,ky))*(1+(ky!=0?kx*kx/ky/ky:0));
+      num+=(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(g[klookup(kx,ky,2*K1+1,2*K2+1)])*(1+(ky!=0?kx*kx/ky/ky:0));
     }
   }
 
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index e23b457..e64eeda 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -13,16 +13,16 @@ typedef struct fft_vects {
 } fft_vect;
 
 // compute u_k
-int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
+int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int nthreads);
 
 // compute the energy as a function of time
-int energy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
+int energy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int nthreads);
 
 // compute enstrophy
-int enstrophy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
+int enstrophy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int nthreads);
 
 // compute solution as a function of time, but do not print anything (useful for debugging)
-int quiet( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double (*g)(int,int), unsigned int nthreads);
+int quiet( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int nthreads);
 
 
 // initialize vectors for computation
@@ -34,10 +34,10 @@ int ns_free_tmps( _Complex double* u, _Complex double* tmp1, _Complex double *tm
 int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2);
 
 // next time step for Irreversible Navier-Stokes equation
-int ins_step( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3);
+int ins_step( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3);
 
 // right side of Irreversible Navier-Stokes equation
-int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2, fft_vect ifft);
+int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft);
 
 // convolution term in right side of equation
 int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft_vect fft2, fft_vect ifft);
@@ -46,7 +46,7 @@ int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft
 int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);
 
 // compute alpha
-_Complex double compute_alpha( _Complex double* u, int K1, int K2, _Complex double (*g)(int,int));
+_Complex double compute_alpha( _Complex double* u, int K1, int K2, _Complex double* g);
 
 // get index for kx,ky in array of size S
 int klookup( int kx, int ky, int S1, int S2);
-- 
cgit v1.2.3-54-g00ecf