From 459ff94855b2659de14ac2ac902b9658d28e19d5 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo@CAE-linux.(none)>
Date: Mon, 21 Apr 2014 11:11:37 +0000
Subject: [PATCH] Added complex spherical harmonics.
---
src/Makefile.am | 10 +-
src/main.c | 3
src/timestep.c | 7 +-
src/tape | 6 +-
src/vesicle.c | 3
src/shcomplex.c | 127 ++++++++++++++++++++++++++++++++++++++++++
src/initial_distribution.c | 3
src/general.h | 2
src/shcomplex.h | 8 ++
9 files changed, 155 insertions(+), 14 deletions(-)
diff --git a/src/Makefile.am b/src/Makefile.am
index 970dec2..7e77543 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,17 +1,17 @@
trisurfdir=../
trisurf_PROGRAMS = trisurf
-trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c bondflip.c main.c poly.c stats.c sh.c
+trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c bondflip.c main.c poly.c stats.c sh.c shcomplex.c
#trisurf_LDFLAGS = -lm -lconfuse
shdiscoverdir=../
shdiscover_PROGRAMS= shdiscover
-shdiscover_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c shdiscover.c poly.c stats.c
+shdiscover_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c shdiscover.c poly.c stats.c shcomplex.c
AM_CFLAGS = -Wall -Werror
co_testdir=../
co_test_PROGRAMS=co_test
-co_test_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c co_test.c frame.c bondflip.c poly.c stats.c
+co_test_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c co_test.c frame.c bondflip.c poly.c stats.c shcomplex.c
spherical_trisurfdir=../
spherical_trisurf_PROGRAMS = spherical_trisurf
-spherical_trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf.c sh.c bondflip.c poly.c stats.c
+spherical_trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf.c sh.c bondflip.c poly.c stats.c shcomplex.c
spherical_trisurf_ffdir=../
spherical_trisurf_ff_PROGRAMS = spherical_trisurf_ff
-spherical_trisurf_ff_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf_ff.c sh.c bondflip.c poly.c stats.c
+spherical_trisurf_ff_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf_ff.c sh.c bondflip.c poly.c stats.c shcomplex.c
diff --git a/src/general.h b/src/general.h
index e47ef1c..8165505 100644
--- a/src/general.h
+++ b/src/general.h
@@ -3,6 +3,7 @@
#include<stdarg.h>
#include<stdio.h>
+#include<gsl/gsl_complex.h>
/* @brief This is a header file, defining general constants and structures.
* @file header.h
* @author Samo Penic
@@ -222,6 +223,7 @@
typedef struct {
ts_uint l;
ts_double **ulm;
+ gsl_complex **ulmComplex;
ts_double **sumUlm2;
ts_uint N;
ts_double **co;
diff --git a/src/initial_distribution.c b/src/initial_distribution.c
index 5006a64..bff9f05 100644
--- a/src/initial_distribution.c
+++ b/src/initial_distribution.c
@@ -12,6 +12,7 @@
#include "poly.h"
#include "io.h"
#include "sh.h"
+#include "shcomplex.h"
ts_vesicle *initial_distribution_dipyramid(ts_uint nshell, ts_uint ncmax1, ts_uint ncmax2, ts_uint ncmax3, ts_double stepsize){
ts_fprintf(stdout,"Starting initial_distribution on vesicle with %u shells!...\n",nshell);
@@ -94,7 +95,7 @@
vesicle->pressure= tape->pressure;
vesicle->pswitch=tape->pswitch;
if(tape->shc>0){
- vesicle->sphHarmonics=sph_init(vesicle->vlist,tape->shc);
+ vesicle->sphHarmonics=complex_sph_init(vesicle->vlist,tape->shc);
}
else {
vesicle->sphHarmonics=NULL;
diff --git a/src/main.c b/src/main.c
index 9bb3430..25f1343 100644
--- a/src/main.c
+++ b/src/main.c
@@ -12,6 +12,7 @@
#include "timestep.h"
#include "poly.h"
#include "sh.h"
+#include "shcomplex.h"
/** Entrance function to the program
* @param argv is a number of parameters used in program call (including the program name
@@ -53,7 +54,7 @@
vesicle->clist->dmin_interspecies = tape->dmin_interspecies*tape->dmin_interspecies;
/* spherical harmonics */
if(tape->shc>0){
- vesicle->sphHarmonics=sph_init(vesicle->vlist,tape->shc);
+ vesicle->sphHarmonics=complex_sph_init(vesicle->vlist,tape->shc);
}
else {
vesicle->sphHarmonics=NULL;
diff --git a/src/shcomplex.c b/src/shcomplex.c
new file mode 100644
index 0000000..d9488ca
--- /dev/null
+++ b/src/shcomplex.c
@@ -0,0 +1,127 @@
+#include<math.h>
+#include<stdlib.h>
+#include<gsl/gsl_complex.h>
+#include<gsl/gsl_complex_math.h>
+#include<gsl/gsl_sf_legendre.h>
+#include "general.h"
+#include "sh.h"
+#include "shcomplex.h"
+
+
+ts_spharm *complex_sph_init(ts_vertex_list *vlist, ts_uint l){
+ ts_uint j,i;
+ ts_spharm *sph=(ts_spharm *)malloc(sizeof(ts_spharm));
+
+ sph->N=0;
+ /* lets initialize Ylm for each vertex. */
+ sph->Ylmi=(ts_double ***)calloc(l,sizeof(ts_double **));
+ for(i=0;i<l;i++){
+ sph->Ylmi[i]=(ts_double **)calloc(2*i+1,sizeof(ts_double *));
+ for(j=0;j<(2*i+1);j++){
+ sph->Ylmi[i][j]=(ts_double *)calloc(vlist->n,sizeof(ts_double));
+ }
+ }
+
+ /* lets initialize ulm */
+ sph->ulm=(ts_double **)calloc(l,sizeof(ts_double *));
+ sph->ulmComplex=(gsl_complex **)calloc(l,sizeof(gsl_complex *));
+ for(j=0;j<l;j++){
+ sph->ulm[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
+ sph->ulmComplex[j]=(gsl_complex *)calloc(2*j+1,sizeof(gsl_complex));
+ }
+
+ /* lets initialize sum of Ulm2 */
+ sph->sumUlm2=(ts_double **)calloc(l,sizeof(ts_double *));
+ for(j=0;j<l;j++){
+ sph->sumUlm2[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
+ }
+
+ /* lets initialize co */
+//NOTE: C is has zero based indexing. Code is imported from fortran and to comply with original indexes we actually generate one index more. Also second dimension is 2*j+2 instead of 2*j+2. elements starting with 0 are useles and should be ignored!
+ sph->co=(ts_double **)calloc(l+1,sizeof(ts_double *));
+ for(j=0;j<=l;j++){
+ sph->co[j]=(ts_double *)calloc(2*j+2,sizeof(ts_double));
+ }
+
+ sph->l=l;
+
+ /* Calculate coefficients that will remain constant during all the simulation */
+ precomputeShCoeff(sph);
+
+ return sph;
+}
+
+ts_bool complex_sph_free(ts_spharm *sph){
+ int i,j;
+ if(sph==NULL) return TS_FAIL;
+ for(i=0;i<sph->l;i++){
+ if(sph->ulm[i]!=NULL) free(sph->ulm[i]);
+ if(sph->ulmComplex[i]!=NULL) free(sph->ulmComplex[i]);
+ if(sph->sumUlm2[i]!=NULL) free(sph->sumUlm2[i]);
+ if(sph->co[i]!=NULL) free(sph->co[i]);
+ }
+ if(sph->co[sph->l]!=NULL) free(sph->co[sph->l]);
+ if(sph->co != NULL) free(sph->co);
+ if(sph->ulm !=NULL) free(sph->ulm);
+ if(sph->ulmComplex !=NULL) free(sph->ulmComplex);
+
+ if(sph->Ylmi!=NULL) {
+ for(i=0;i<sph->l;i++){
+ if(sph->Ylmi[i]!=NULL){
+ for(j=0;j<i*2+1;j++){
+ if(sph->Ylmi[i][j]!=NULL) free (sph->Ylmi[i][j]);
+ }
+ free(sph->Ylmi[i]);
+ }
+ }
+ free(sph->Ylmi);
+ }
+
+ free(sph);
+ return TS_SUCCESS;
+}
+
+
+ts_bool calculateUlmComplex(ts_vesicle *vesicle){
+ ts_int i,j,k,m,l;
+ ts_vertex *cvtx;
+ ts_coord coord;
+/* set all values to zero */
+ for(i=0;i<vesicle->sphHarmonics->l;i++){
+ for(j=0;j<2*i+1;j++) GSL_SET_COMPLEX(&(vesicle->sphHarmonics->ulmComplex[i][j]),0.0,0.0);
+ }
+
+ for(k=0;k<vesicle->vlist->n; k++){
+ cvtx=vesicle->vlist->vtx[k];
+ cart2sph(&coord,cvtx->x,cvtx->y,cvtx->z);
+ for(i=0;i<vesicle->sphHarmonics->l;i++){
+ for(j=0;j<2*i+1;j++){
+ m=j-i;
+ l=i;
+ if(m>=0){
+ // fprintf(stderr, "Racunam za l=%d, m=%d\n", l,m);
+ vesicle->sphHarmonics->ulmComplex[i][j]=gsl_complex_add(vesicle->sphHarmonics->ulmComplex[i][j], gsl_complex_conjugate(gsl_complex_mul_real(gsl_complex_polar(1.0,(ts_double)m*coord.e2),cvtx->solAngle*cvtx->relR*gsl_sf_legendre_sphPlm(l,m,cos(coord.e3)))) );
+ } else {
+ // fprintf(stderr, "Racunam za l=%d, abs(m=%d)\n", l,m);
+ vesicle->sphHarmonics->ulmComplex[i][j]=gsl_complex_add(vesicle->sphHarmonics->ulmComplex[i][j], gsl_complex_conjugate(gsl_complex_mul_real(gsl_complex_polar(1.0,(ts_double)m*coord.e2),cvtx->solAngle*cvtx->relR*pow(-1,m)*gsl_sf_legendre_sphPlm(l,-m,cos(coord.e3)))) );
+
+ }
+ }
+ }
+ }
+ return TS_SUCCESS;
+}
+
+ts_bool storeUlmComplex2(ts_vesicle *vesicle){
+
+ ts_spharm *sph=vesicle->sphHarmonics;
+ ts_int i,j;
+ for(i=0;i<sph->l;i++){
+ for(j=0;j<2*i+1;j++){
+ sph->sumUlm2[i][j]+=gsl_complex_abs2(sph->ulmComplex[i][j]);
+ }
+ }
+ sph->N++;
+ return TS_SUCCESS;
+}
+
diff --git a/src/shcomplex.h b/src/shcomplex.h
new file mode 100644
index 0000000..161be6a
--- /dev/null
+++ b/src/shcomplex.h
@@ -0,0 +1,8 @@
+#ifndef _H_SH_COMPLEX
+#define _H_SH_COMPLEX
+#include "general.h"
+ts_bool storeUlmComplex2(ts_vesicle *vesicle);
+ts_spharm *complex_sph_init(ts_vertex_list *vlist, ts_uint l);
+ts_bool complex_sph_free(ts_spharm *sph);
+ts_bool calculateUlmComplex(ts_vesicle *vesicle);
+#endif
diff --git a/src/tape b/src/tape
index 7de3646..be15cfd 100644
--- a/src/tape
+++ b/src/tape
@@ -44,15 +44,15 @@
####### Program Control ############
#how many MC sweeps between subsequent records of states to disk
-mcsweeps=10
+mcsweeps=250
#how many initial mcsweeps*inititer MC sweeps before recording to disk?
inititer=0
#how many records do you want on the disk iteration are there in a run?
-iterations=10
+iterations=1000
###### Spherical harmonics ###########
-spherical_harmonics_coefficients=21
+spherical_harmonics_coefficients=15
#shut up if we are using cluster!!!
quiet=false
diff --git a/src/timestep.c b/src/timestep.c
index f24eccd..62e640c 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -10,6 +10,7 @@
#include "io.h"
#include "stats.h"
#include "sh.h"
+#include "shcomplex.h"
#include "vesicle.h"
ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations, ts_uint start_iteration){
@@ -51,9 +52,9 @@
r0=getR0(vesicle);
if(vesicle->sphHarmonics!=NULL){
preparationSh(vesicle,r0);
- calculateYlmi(vesicle);
- calculateUlm(vesicle);
- storeUlm2(vesicle);
+ //calculateYlmi(vesicle);
+ calculateUlmComplex(vesicle);
+ storeUlmComplex2(vesicle);
saveAvgUlm2(vesicle);
fd1=fopen("state.dat","w");
fprintf(fd1,"%e %e\n",vesicle->volume, getR0(vesicle));
diff --git a/src/vesicle.c b/src/vesicle.c
index 0a1b84b..f63eb12 100644
--- a/src/vesicle.c
+++ b/src/vesicle.c
@@ -7,6 +7,7 @@
#include "stdlib.h"
#include "poly.h"
#include "sh.h"
+#include "shcomplex.h"
ts_vesicle *init_vesicle(ts_uint N, ts_uint ncmax1, ts_uint ncmax2, ts_uint
ncmax3, ts_double stepsize){
@@ -36,7 +37,7 @@
triangle_list_free(vesicle->tlist);
cell_list_free(vesicle->clist);
poly_list_free(vesicle->poly_list);
- sph_free(vesicle->sphHarmonics);
+ complex_sph_free(vesicle->sphHarmonics);
free(vesicle);
return TS_SUCCESS;
}
--
Gitblit v1.9.3