Commit 618b5766 authored by Olivier Bonnefon's avatar Olivier Bonnefon
Browse files

add adaptive ts simulator

parent c86fb782
#include <iostream>
#include <cfloat>
using namespace std;
#include "error.hpp"
#include "AFunction.hpp"
#include "rgraph.hpp"
#include "RNM.hpp"
#include "fem.hpp"
#include "FESpace.hpp"
#include "MeshPoint.hpp"
#include "AFunction_ext.hpp" // Extension of "AFunction.hpp" to deal with more than 3 parameters function
using namespace Fem2D;
double SaveVecAppend(KN<double> *const & f, string *const & nome){
std::ofstream outfile (nome->data(),ios_base::binary | std::ofstream::app);
long int nn = f->N(); // get number of nodes
long int dim=nn;
long pos = outfile.tellp();
if (!pos)
outfile.write ((char*) &dim, sizeof(long int));//write the dimension of the vector
double ftemp ;
for(long int i=0; i<nn; i++) {
ftemp = *(f[0]+i) ;
outfile.write ((char*) &ftemp, sizeof(double));
}
outfile.close();
return 0.0; // dummy return value.
}
double SetSizeVec(long int *const & size, string *const & nome){
//std::ofstream outfile (nome->data(),ios_base::binary | std::ofstream::app);
std::fstream outfile(nome->data(), ios_base::binary | std::fstream::out | std::fstream::in);
long int N=*size;
outfile.write ((char*) &N, sizeof(long int));//write the dimension of the vector
outfile.close();
return 0.0; // dummy return value.
}
double SaveVec(KN<double> *const & f, string *const & nome)
{
std::ofstream outfile (nome->data(),ios_base::binary);
// To access value at node i of vector N, do as follow: *(N[0]+i)
// Explanation (C++ for dummies as I am ;-):
// N is an alias to the KN object.
// N[0] is a pointer to the first element of the vector.
// N[0]+i is a pointer to the ith element of the vector.
// *(N[0]+i) is the value of the ith element of the vector.
long int nn = f->N(); // get number of nodes
long int dim=nn;
outfile.write ((char*) &dim, sizeof(long int));//write the dimension of the vector
double ftemp ;
for(long int i=0; i<nn; i++) {
ftemp = *(f[0]+i) ;
outfile.write ((char*) &ftemp, sizeof(double));
}
outfile.close();
return 0.0; // dummy return value.
}
double GetSizeVec(long int *const & size,string *const & nome){
std::ifstream infile (nome->data(),ios_base::binary);
long int dim;
infile.read((char *) &dim, sizeof(long int));
*size=dim;
infile.close();
return 0.0;
}
double LoadVec(KN<double> *const & ww, string *const & nome)
{
std::ifstream infile (nome->data(),ios_base::binary);
long int dim;
infile.read((char *) &dim, sizeof(long int));
double dtemp;
for(long int i=0; i<dim; i++)
{
infile.read((char *) &dtemp, sizeof(double));
*(ww[0]+i)=dtemp ;
}
return 0.0; // dummy return value.
}
double LoadFlag(long int *const & ww, string *const & nome)
{
std::ifstream infile (nome->data(),ios_base::binary);
long int flag;
infile.read((char *) &flag, sizeof(long int));
*ww=flag;
return 0.0; // dummy return value.
}
double flag(long int *const & FLAG,string *const &nome)
{
std::ofstream outfile (nome->data(),ios_base::binary);
long int Flag;
Flag= *FLAG;
outfile.write ((char*) &Flag, sizeof(long int));
outfile.close();
return 0.0;
}
// add the function name to the freefem++ table
/* class Init { public:
Init();
};
$1 */
static void Load_Init(){
Global.Add("LoadVec","(",new OneOperator2_<double, KN<double>*, string* >(LoadVec));
Global.Add("LoadFlag","(",new OneOperator2_<double,long int*, string* >(LoadFlag));
Global.Add("SaveVec","(",new OneOperator2_<double,KN<double>*, string* >(SaveVec));
Global.Add("flag","(",new OneOperator2_<double,long int*,string* >(flag));
Global.Add("GetSizeVec","(",new OneOperator2_<double,long int*,string* >(GetSizeVec));
Global.Add("SaveVecAppend","(",new OneOperator2_<double,KN<double>*, string* >(SaveVecAppend));
Global.Add("SetSizeVec","(",new OneOperator2_<double,long int*, string* >(SetSizeVec));
}
LOADFUNC(Load_Init)
load "BinaryIO"
// test of BinaryIO tools ...
real[int] buffer1(5);
real[int] buffer2(5);
for (int ii=0;ii<5;ii++)
buffer1[ii]=ii;
SaveVec(buffer1,"saveBuff");
int size=0;
GetSizeVec(size,"saveBuff");
cout<<"size ="<<size<<endl;
LoadVec(buffer2,"saveBuff");
for (int ii=0;ii<size;ii++)
cout<<buffer2[ii]<<" ";
cout<<endl;
for (int jj=0;jj<4;jj++){
for (int ii=0;ii<5;ii++)
buffer1[ii]=buffer1[ii]+5;
SaveVecAppend(buffer1,"saveBuff");
}
int newSize=25;
SetSizeVec(newSize,"saveBuff");
GetSizeVec(size,"saveBuff");
real[int] buffer3(size);
cout<<"now size ="<<size<<endl;
LoadVec(buffer3,"saveBuff");
for (int ii=0;ii<size;ii++)
cout<<buffer3[ii]<<" ";
cout<<endl;
load "BinaryIO"
int size;
GetSizeVec(size,"SAVE/u2d0__rep_0");
cout<<"now size ="<<size<<endl;
This file is copied during the execution of the C prog gen2speciesNL.c.
It contains the command that must be done to allow the simulation.
This commands are necessary to run the simulation:
1) more userFDef0.txt
Just to check the dynamic on each patch
2) chmod +x genSetFunc.sh
3) ./genSetFunc.sh
To generate the dynamic on each patch
The simuated must be ready to run:
FreeFem++ cas1.edp
Finaly, visualized with 'python vtktopng.py'
Surf0=1.36556e+11
Surf1=3.34226e+10
Surf2=5.28336e+10
Surf3=2.35296e+10
Surf4=4.77306e+10
Surf5=4.20849e+10
Surf6=3.22273e+10
Surf7=1.99198e+10
Surf8=4.59006e+10
Surf9=2.2696e+10
Surf10=2.2557e+10
Surf11=2.88487e+10
Surf12=5.73733e+09
Surf13=3.18954e+09
Surf14=2.29261e+10
ce repeertoire contient les scripts de post tratement de simulation.
1) checkSimu.edp : verifie que la presence et la taille des fichiers de sortie des 100 iterations du simulateur.
2) postT.cpp : calcul les trajectoire min et max
3) visuSimuMax.edp: creer les sorties VTK des trajectoires min max, doit etre lancer une fois avec "max" et une fois avec "min"
4)vtktopng.py et vtktopngGen.py: creer les images a partir des sorties vtk.
Attention, il convient de s'assurer que les scrits et programme pointe bien vers le repertoire de sorties que simulations visé.
\ No newline at end of file
load "BinaryIO";
include "paramsTS.edp";
//string BASEREP="../SAVE/MEDIAN/";
string BASEREP="../SAVE/HIGHT/";
int curStep=0;
verbosity=0;
int[int] sizeBinary(15);
for (int ii=0;ii<15;ii++)
sizeBinary[ii]=0;
int sizeBinaryAux=0;
macro checkFile(XXINT){
string FILENAME=BASEREP+"u2d"+string(XXINT)+"__rep_"+string(numRun);
if (exec("ls "+FILENAME+" 2>/dev/null >/dev/null")){
nofile=1;
}else{
if (!sizeBinary[XXINT]){
GetSizeVec(sizeBinary[XXINT],FILENAME);
}else{
GetSizeVec(sizeBinaryAux,FILENAME);
if (sizeBinaryAux != sizeBinary[XXINT])
cout<<"ERROR "+FILENAME+" size is "<<sizeBinaryAux<<" and should be "<<sizeBinary[XXINT]<<endl;
}
}
}//
for( int numRun=0;numRun<100;numRun++){
int nofile=0;
//for (CURYEAR=FIRSTYEAR; CURYEAR<=LASTYEAR; CURYEAR++){
// for ( curStep=0;curStep<nbSteps;curStep++){
for (int ii=0;ii<15;ii++)
checkFile(ii);
if (!nofile){
cout<<"numRun="<<numRun<<" ok"<<endl;
}else{
cout<<"numRun="<<numRun<<" missing"<<endl;
}
}
includepath += ".."
includepath += "../../Tools1D/"
includepath += "../../Recorder/"
includepath += "../../examples1D/"
loadpath += "/opt/freefem++-3.58/examples++-load/"
#include <fstream>
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
using namespace std;
#define BASEREP "../SAVE/HIGHT/"
long int GetSizeVec(char* name){
std::ifstream infile (name,ios_base::binary);
long int dim;
infile.read((char *) &dim, sizeof(long int));
infile.close();
return dim;
}
void LoadVec(char* name,double *pD)
{
std::ifstream infile (name,ios_base::binary);
long int dim;
infile.read((char *) &dim, sizeof(long int));
double dtemp;
for(long int i=0; i<dim; i++)
{
infile.read((char *) &dtemp, sizeof(double));
*(pD++)=dtemp ;
}
infile.close();
}
void SaveVec(double * pD, long int N, char* name)
{
std::ofstream outfile (name,ios_base::binary);
outfile.write ((char*) &N, sizeof(long int));//write the dimension of the vector
double ftemp ;
for(long int i=0; i<N; i++) {
ftemp = pD[i] ;
outfile.write ((char*) &ftemp, sizeof(double));
}
outfile.close();
}
int main(){
char filename[256];
int N=4;
double * pBufTraj[15];
double * pMaxTraj[15];
double * pMinTraj[15];
long int size[15];
for (int j=0;j<15;j++){
sprintf(filename,"%su2d%i__rep_0",BASEREP,j);
size[j]=GetSizeVec(filename);
pBufTraj[j]=(double *)malloc(size[j]*sizeof(double));
pMaxTraj[j]=(double *)malloc(size[j]*sizeof(double));
pMinTraj[j]=(double *)malloc(size[j]*sizeof(double));
LoadVec(filename,pMaxTraj[j]);
memcpy(pMinTraj[j],pMaxTraj[j],size[j]*sizeof(double));
}
for(int numRun=1;numRun<N;numRun++){
for (int j=0;j<15;j++){
sprintf(filename,"%su2d%i__rep_%i",BASEREP,j,numRun);
LoadVec(filename,pBufTraj[j]);
for (int k=0;k<size[j];k++){
if (pBufTraj[j][k]>pMaxTraj[j][k])
pMaxTraj[j][k]=pBufTraj[j][k];
if (pBufTraj[j][k]<pMinTraj[j][k])
pMinTraj[j][k]=pBufTraj[j][k];
}
}
}
for (int j=0;j<15;j++){
sprintf(filename,"%su2d%i__max",BASEREP,j);
SaveVec(pMaxTraj[j],size[j],filename);
sprintf(filename,"%su2d%i__min",BASEREP,j);
SaveVec(pMinTraj[j],size[j],filename);
}
for (int j=0;j<15;j++){
free(pBufTraj[j]);
free(pMaxTraj[j]);
free(pMinTraj[j]);
}
return 0;
}
savevtk(OUTPUTDIR+"D2/u"+maxormin+"0s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th0,utm12d0,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"1s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th1,utm12d1,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"2s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th2,utm12d2,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"3s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th3,utm12d3,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"4s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th4,utm12d4,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"5s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th5,utm12d5,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"6s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th6,utm12d6,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"7s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th7,utm12d7,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"8s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th8,utm12d8,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"9s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th9,utm12d9,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"10s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th10,utm12d10,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"11s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th11,utm12d11,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"12s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th12,utm12d12,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"13s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th13,utm12d13,dataname="UU");
savevtk(OUTPUTDIR+"D2/u"+maxormin+"14s"+string(CURYEAR)+"_"+string(curStep)+".vtk",Th14,utm12d14,dataname="UU");
This diff is collapsed.
load "BinaryIO";
include "paramsTS.edp";
include "myMacro.edp";
include "defGeomGen.edp"
include "defDofsGen.edp"
mesh ThALL=readmesh("../france.msh");
plot(ThALL,wait=true);
fespace VALL(ThALL,P1);
VALL vall;
VALL vabuf;
matrix IVA0=interpolate(VALL,Vh0,inside=1);
matrix IVA1=interpolate(VALL,Vh1,inside=1);
matrix IVA2=interpolate(VALL,Vh2,inside=1);
matrix IVA3=interpolate(VALL,Vh3,inside=1);
matrix IVA4=interpolate(VALL,Vh4,inside=1);
matrix IVA5=interpolate(VALL,Vh5,inside=1);
matrix IVA6=interpolate(VALL,Vh6,inside=1);
matrix IVA7=interpolate(VALL,Vh7,inside=1);
matrix IVA8=interpolate(VALL,Vh8,inside=1);
matrix IVA9=interpolate(VALL,Vh9,inside=1);
matrix IVA10=interpolate(VALL,Vh10,inside=1);
matrix IVA11=interpolate(VALL,Vh11,inside=1);
matrix IVA12=interpolate(VALL,Vh12,inside=1);
matrix IVA13=interpolate(VALL,Vh13,inside=1);
matrix IVA14=interpolate(VALL,Vh14,inside=1);
int curStep=0;
int numRun=0;
int size;
//mesh ThALL=readmesh("../france.msh");
//fespace VALL(ThALL,P1);
fespace VALL0(ThALL,P0);
VALL UALL;
UALL=0;
plot(ThALL,UALL);
VALL0 nuT;
for (int i=0;i<ThALL.nt;i++)
nuT[][i]=i;
int NlocationsLamberts=2;
real[int] Xlambert(NlocationsLamberts);
real[int] Ylambert(NlocationsLamberts);
int[int] numDoflambert(NlocationsLamberts);
//coordonne avignon en lambert II etendu
Xlambert[0]=798256;
Ylambert[0]=1886174;
//coordonne paris en lambert II etendu
Xlambert[1]=601638;
Ylambert[1]=2428265;
for (int i=0;i<NlocationsLamberts;i++){
int numT=nuT(Xlambert[i],Ylambert[i]);
numDoflambert[i]=ThALL[numT][0];
UALL[][numDoflambert[i]]=1;
}
plot(ThALL,UALL);
string REPSIMU="../SAVE/LOW/";
ofstream pTraj("traj.txt");
{
GetSizeVec(size,REPSIMU+"u2d0__rep_"+string(numRun));
real[int] u0(size);
LoadVec(u0,REPSIMU+"u2d0__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d1__rep_"+string(numRun));
real[int] u1(size);
LoadVec(u1,REPSIMU+"u2d1__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d2__rep_"+string(numRun));
real[int] u2(size);
LoadVec(u2,REPSIMU+"u2d2__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d3__rep_"+string(numRun));
real[int] u3(size);
LoadVec(u3,REPSIMU+"u2d3__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d4__rep_"+string(numRun));
real[int] u4(size);
LoadVec(u4,REPSIMU+"u2d4__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d5__rep_"+string(numRun));
real[int] u5(size);
LoadVec(u5,REPSIMU+"u2d5__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d6__rep_"+string(numRun));
real[int] u6(size);
LoadVec(u6,REPSIMU+"u2d6__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d7__rep_"+string(numRun));
real[int] u7(size);
LoadVec(u7,REPSIMU+"u2d7__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d8__rep_"+string(numRun));
real[int] u8(size);
LoadVec(u8,REPSIMU+"u2d8__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d9__rep_"+string(numRun));
real[int] u9(size);
LoadVec(u9,REPSIMU+"u2d9__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d10__rep_"+string(numRun));
real[int] u10(size);
LoadVec(u10,REPSIMU+"u2d10__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d11__rep_"+string(numRun));
real[int] u11(size);
LoadVec(u11,REPSIMU+"u2d11__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d12__rep_"+string(numRun));
real[int] u12(size);
LoadVec(u12,REPSIMU+"u2d12__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d13__rep_"+string(numRun));
real[int] u13(size);
LoadVec(u13,REPSIMU+"u2d13__rep_"+string(numRun));
GetSizeVec(size,REPSIMU+"u2d14__rep_"+string(numRun));
real[int] u14(size);
LoadVec(u14,REPSIMU+"u2d14__rep_"+string(numRun));
int nbCallSave=0;
for (CURYEAR=FIRSTYEAR; CURYEAR<=LASTYEAR; CURYEAR++){
int firstStep=0;
if (CURYEAR==2011)
firstStep=362-122;
for ( curStep=firstStep;curStep<nbSteps-1;curStep++){
cout<<CURYEAR<<" "<<curStep<<endl;
utm12d0[]=u0(nbCallSave*utm12d0.n:(nbCallSave+1)*utm12d0.n);
utm12d1[]=u1(nbCallSave*utm12d1.n:(nbCallSave+1)*utm12d1.n);
utm12d2[]=u2(nbCallSave*utm12d2.n:(nbCallSave+1)*utm12d2.n);
utm12d3[]=u3(nbCallSave*utm12d3.n:(nbCallSave+1)*utm12d3.n);
utm12d4[]=u4(nbCallSave*utm12d4.n:(nbCallSave+1)*utm12d4.n);
utm12d5[]=u5(nbCallSave*utm12d5.n:(nbCallSave+1)*utm12d5.n);
utm12d6[]=u6(nbCallSave*utm12d6.n:(nbCallSave+1)*utm12d6.n);
utm12d7[]=u7(nbCallSave*utm12d7.n:(nbCallSave+1)*utm12d7.n);
utm12d8[]=u8(nbCallSave*utm12d8.n:(nbCallSave+1)*utm12d8.n);
utm12d9[]=u9(nbCallSave*utm12d9.n:(nbCallSave+1)*utm12d9.n);
utm12d10[]=u10(nbCallSave*utm12d10.n:(nbCallSave+1)*utm12d10.n);
utm12d11[]=u11(nbCallSave*utm12d11.n:(nbCallSave+1)*utm12d11.n);
utm12d12[]=u12(nbCallSave*utm12d12.n:(nbCallSave+1)*utm12d12.n);
utm12d13[]=u13(nbCallSave*utm12d13.n:(nbCallSave+1)*utm12d13.n);
utm12d14[]=u14(nbCallSave*utm12d14.n:(nbCallSave+1)*utm12d14.n);
VALL vabuf=0;
VALL vall=0;
vabuf[]=IVA0*utm12d0[];
adddofs(vall,vabuf);
vabuf[]=IVA1*utm12d1[];
adddofs(vall,vabuf);
vabuf[]=IVA2*utm12d2[];
adddofs(vall,vabuf);
vabuf[]=IVA3*utm12d3[];
adddofs(vall,vabuf);
vabuf[]=IVA4*utm12d4[];
adddofs(vall,vabuf);
vabuf[]=IVA5*utm12d5[];
adddofs(vall,vabuf);
vabuf[]=IVA6*utm12d6[];
adddofs(vall,vabuf);
vabuf[]=IVA7*utm12d7[];
adddofs(vall,vabuf);
vabuf[]=IVA8*utm12d8[];
adddofs(vall,vabuf);
vabuf[]=IVA9*utm12d9[];
adddofs(vall,vabuf);
vabuf[]=IVA10*utm12d10[];
adddofs(vall,vabuf);
vabuf[]=IVA11*utm12d11[];
adddofs(vall,vabuf);
vabuf[]=IVA12*utm12d12[];
adddofs(vall,vabuf);
vabuf[]=IVA13*utm12d13[];
adddofs(vall,vabuf);
vabuf[]=IVA14*utm12d14[];
adddofs(vall,vabuf);
for (int i=0;i<NlocationsLamberts;i++){
pTraj<<vall(Xlambert[i],Ylambert[i])<<" ";
}
pTraj<<"\n";
if (0 && !(nbCallSave % 10)){
plot(utm12d0,utm12d1,utm12d2,utm12d3,utm12d4,utm12d5,utm12d6,utm12d7,utm12d8,utm12d9,utm12d10,utm12d11,utm12d12,utm12d13,utm12d14,
cmm="rep "+string(numRun)+" year="+string(CURYEAR)+" step="+string(curStep),value=true,fill=true,wait=false);
plot(vall,cmm="all",value=true,fill=true,wait=false);
}
nbCallSave++;
}