From 28a4a3879c8e5ca8143ec1a7f845394fa5e1674d Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Tue, 28 May 2024 16:44:50 +0200
Subject: [PATCH 01/23] add a new class, wich is a subclass of model and is
used for second order term with a txt file
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 39 +++++++++++++++++++++
1 file changed, 39 insertions(+)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index ae11721..b0cdb7e 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -42,3 +42,42 @@ class diffusionModel(model):
def updateBlockMatAdjoint(self,compi,compj):
if (compi==compj):
return self.blockAdjoint
+
+class diffusionModelTxt(model):
+ def __init__(self,aDs, aTxtFile) -> None:
+ super().__init__(aDs)
+ if (exists(aTxtFile+".txt")):
+ self.txtFile=aTxtFile
+ else:
+ self.txtFile=computerEnv.modelsDir+aTxtFile #aTxtFile form is "affineGrowth"
+ model.myPrint("Looking for model "+self.txtFile)
+ try:
+ f = open(self.txtFile+".txt", "r")
+ except IOError:
+ model.myPrint("IOError in txtModel when reading "+self.txtFile+".txt")
+ return
+ try:
+ fd = open(self.txtFile+"DVP.txt", "r")
+ except FileNotFoundError:
+ self.dvFileExist=False
+ print("Création fichier "+self.txtFile+"DVP.tx")
+ fd = open(self.txtFile+"Dv.txt", "w") #open file in mod write
+ for numComps in range(self.nbComps):
+ strF=f.readline()
+ strF=strF.split('=')[1].strip()#.replace("^","**")
+ for p in range(self.aDs.study.nbParam):
+ grad1s=diff(strF,'p'+str(p+1))
+ grad1=str(grad1s).replace("**","^")
+ #write in fic
+ fd.write("d"+txtModel.unknownName[k]+"f"+txtModel.unknownName[numComps]+"="+grad1+"\n")
+ f.close()
+ fd.close()
+ f = open(self.txtFile+".txt", "r")
+ try:
+ fd = open(self.txtFile+"Dv.txt", "r")
+ except IOError:
+ model.myPrint("IOError in txtModel when reading "+self.txtFile+"Dv.txt")
+ return
+ except IOError:
+ model.myPrint("IOError in txtModel when reading "+self.txtFile+"Dv.txt")
+ return
--
GitLab
From f1d37bd36341fc48dea175fcfe26195092aaf1e6 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Tue, 28 May 2024 17:05:50 +0200
Subject: [PATCH 02/23] verifies existence of file and create derivate file if
did'nt exist
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 6 ++++++
1 file changed, 6 insertions(+)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index b0cdb7e..7149d70 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -44,8 +44,12 @@ class diffusionModel(model):
return self.blockAdjoint
class diffusionModelTxt(model):
+ """ Class used to create a second order term with a txt file
+
+ """
def __init__(self,aDs, aTxtFile) -> None:
super().__init__(aDs)
+ # we check if file exist
if (exists(aTxtFile+".txt")):
self.txtFile=aTxtFile
else:
@@ -56,9 +60,11 @@ class diffusionModelTxt(model):
except IOError:
model.myPrint("IOError in txtModel when reading "+self.txtFile+".txt")
return
+ # we check if the derivate of parameter file exist
try:
fd = open(self.txtFile+"DVP.txt", "r")
except FileNotFoundError:
+ # if didn't exist, we create it
self.dvFileExist=False
print("Création fichier "+self.txtFile+"DVP.tx")
fd = open(self.txtFile+"Dv.txt", "w") #open file in mod write
--
GitLab
From 32181870edb53feb8141eeffa11938c768ffbbf8 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Tue, 28 May 2024 17:46:05 +0200
Subject: [PATCH 03/23] patch name of file
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 11 ++++++-----
1 file changed, 6 insertions(+), 5 deletions(-)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 7149d70..3a3f77e 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -67,7 +67,7 @@ class diffusionModelTxt(model):
# if didn't exist, we create it
self.dvFileExist=False
print("Création fichier "+self.txtFile+"DVP.tx")
- fd = open(self.txtFile+"Dv.txt", "w") #open file in mod write
+ fd = open(self.txtFile+"DVP.txt", "w") #open file in mod write
for numComps in range(self.nbComps):
strF=f.readline()
strF=strF.split('=')[1].strip()#.replace("^","**")
@@ -75,15 +75,16 @@ class diffusionModelTxt(model):
grad1s=diff(strF,'p'+str(p+1))
grad1=str(grad1s).replace("**","^")
#write in fic
- fd.write("d"+txtModel.unknownName[k]+"f"+txtModel.unknownName[numComps]+"="+grad1+"\n")
+ fd.write("dp"+str(p+1)+"f"+txtModel.unknownName[numComps]+"="+grad1+"\n")
f.close()
fd.close()
f = open(self.txtFile+".txt", "r")
try:
- fd = open(self.txtFile+"Dv.txt", "r")
+ fd = open(self.txtFile+"DVP.txt", "r")
except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+"Dv.txt")
+ model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
return
except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+"Dv.txt")
+ model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
return
+ ### find idea for derivate in obprocess.computeGradV
--
GitLab
From 3f748d93e44047912b8b31efa4aaa47d25b17622 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Thu, 30 May 2024 10:25:09 +0200
Subject: [PATCH 04/23] quick update, not finish
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 29 ++++++++++++++-------
1 file changed, 20 insertions(+), 9 deletions(-)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 3a3f77e..858b8e7 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -11,14 +11,8 @@ class diffusionModel(model):
super().__init__(aDs)
### NOTE Test: we define diffusionModel with param
aDs.addSecondOrderTerm(self)
- self.coefName=aDiffCoef
- if type(aDiffCoef) != str:
- ### we supposed that aDiffCoef is like "p?" with ?=int (p1,p2,p3...)
- #NOTE: parameter = dictionnaire ????!!!!
- self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
- else:
- numParam=int(aDiffCoef[1])-1
- self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(self.dS.study.param[numParam]))
+ self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
+ self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
self.block=None
self.blockAdjoint=None
self.bilinearForm=None
@@ -43,7 +37,7 @@ class diffusionModel(model):
if (compi==compj):
return self.blockAdjoint
-class diffusionModelTxt(model):
+class txtDiffusionModel(model):
""" Class used to create a second order term with a txt file
"""
@@ -88,3 +82,20 @@ class diffusionModelTxt(model):
model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
return
### find idea for derivate in obprocess.computeGradV
+ ### TODO: read file + create varfs like txtmodel
+ def truc():
+ pass
+ def _computeBlock(self):
+ aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
+ self.bilinearForm=form(aux)
+ self.block=matrix(assemble_matrix(self.bilinearForm))
+ def updateBlockMat(self,compi,compj):
+ if (compi==compj):
+ return self.block
+ def _computeBlockAdjoint(self):
+ aux=dot(grad(self.dS.u),grad(self.diffCoef*self.dS.v))*self.dS.dx
+ self.bilinearForm=form(aux)
+ self.blockAdjoint=matrix(assemble_matrix(self.bilinearForm))
+ def updateBlockMatAdjoint(self,compi,compj):
+ if (compi==compj):
+ return self.blockAdjoint
--
GitLab
From 41516e8ba7bac15fafa1a91ae3973aba5775adc5 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Thu, 30 May 2024 17:27:45 +0200
Subject: [PATCH 05/23] modif updateBlockMat + add exec compile like txtModel
(to finish)
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 46 +++++++++++++++++++--
1 file changed, 42 insertions(+), 4 deletions(-)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 858b8e7..db1f753 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -83,15 +83,53 @@ class txtDiffusionModel(model):
return
### find idea for derivate in obprocess.computeGradV
### TODO: read file + create varfs like txtmodel
- def truc():
- pass
+ ###Â we have \Delta H(u),
+ ### H_comps, Dp_{ai}H_comps
+ self.varfs=[None]*(self.nbComps*(self.dS.study.nbParam+1))
+ self.varfsTm1=[None]*self.nbComps
+ #for the comp number numComp \in {0,..,self.nbComps-1}:
+ #self.varfs[numComp*(self.nbComps+1)] contains the varf F_numComp for the rhs
+ #self.varfs[numComp*(self.nbComps+1)+(j+1)] contains the derivative d_jF_numComp, j \in {0,..,self.nbComps-1}
+ strkeys=r"[\s\+\-\*/\=\(\)\^]+"
+ for numComps in range(self.nbComps):
+ ##READ SECOND ORDER TERM
+ strline=f.readline().split('=')[1].strip()
+ strline=strForDolf(strline, self.unknownName)
+ strline=strline.replace('^','**')
+ ##replace param
+ for count_p in range(self.dS.paramD):
+ strline=strline.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
+ model.myPrint(" read txt model F"+txtModel.unknownName[numComps]+"="+strline)
+ #### replace COVARIABLES
+ strline=updadeVarfWithCovariables(strline,self.dS.covariables)
+ model.myPrint(strline)
+ #### BUILD VARF USING FK
+ code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+"*self.dS.v*self.dS.dx)"
+ exec(compile(code, 'sumstring', 'exec'))
+ #### BUILD VARF USING UTM1
+ code="self.varfsTm1["+str(numComps)+"]="+self._subFXForm(strline,"utm1")+"*self.dS.v*self.dS.dx)"
+ model.myPrint(code)
+ exec(compile(code, 'sumstring', 'exec'))
+ ## READ DERIVATIVES PARAM
+ #TODO
+ for countP in range(self.dS.paramD):
+ strline=fd.readline().split('=')[1].strip()
+ strline=strForDolf(strline,self.unknownName)
+ #strline=strline.replace('^','**')
+ ##replace param
+ for tmpCountP in enumerate(self.dS.paramD):
+ strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountp)+']')
+ strline=updadeVarfWithCovariables(strline,self.dS.covariables)
+ model.myPrint("txtModel: read txt model d"+txtModel.unknownName[k]+"F"+txtModel.unknownName[numComps]+"="+strline)
+ code="self.varfs["+str(numComps)+"*(self.dS.study.nbParam+1)+"+str(k)+"+1]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
+ exec(compile(code, 'sumstring', 'exec'))
+
def _computeBlock(self):
aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
self.bilinearForm=form(aux)
self.block=matrix(assemble_matrix(self.bilinearForm))
def updateBlockMat(self,compi,compj):
- if (compi==compj):
- return self.block
+ return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)+compj+1]))
def _computeBlockAdjoint(self):
aux=dot(grad(self.dS.u),grad(self.diffCoef*self.dS.v))*self.dS.dx
self.bilinearForm=form(aux)
--
GitLab
From dcf7fec702ca288286988f1fb1027b798013d23a Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Fri, 31 May 2024 10:54:34 +0200
Subject: [PATCH 06/23] add import function strForDolf form file model.py
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 4 +---
1 file changed, 1 insertion(+), 3 deletions(-)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index db1f753..60ec276 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -1,4 +1,4 @@
-from mse.MODELS.models import model
+from mse.MODELS.models import model, strForDolf
import mse.dynamicalSystem
from dolfinx.fem import Constant,form
from petsc4py.PETSc import ScalarType
@@ -111,7 +111,6 @@ class txtDiffusionModel(model):
model.myPrint(code)
exec(compile(code, 'sumstring', 'exec'))
## READ DERIVATIVES PARAM
- #TODO
for countP in range(self.dS.paramD):
strline=fd.readline().split('=')[1].strip()
strline=strForDolf(strline,self.unknownName)
@@ -123,7 +122,6 @@ class txtDiffusionModel(model):
model.myPrint("txtModel: read txt model d"+txtModel.unknownName[k]+"F"+txtModel.unknownName[numComps]+"="+strline)
code="self.varfs["+str(numComps)+"*(self.dS.study.nbParam+1)+"+str(k)+"+1]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
exec(compile(code, 'sumstring', 'exec'))
-
def _computeBlock(self):
aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
self.bilinearForm=form(aux)
--
GitLab
From ce80fb0208ffdef448c4493c16c42f5479f05e36 Mon Sep 17 00:00:00 2001
From: Olivier Bonnefon <olivier.bonnefon@inra.fr>
Date: Wed, 29 May 2024 15:39:41 +0200
Subject: [PATCH 07/23] add covariable in varf, not tested during simulation
---
.../mseMatTools/Frontend/test_MseMatTools.py | 11 ++-
.../Frontend/test_MseMatToolsRaster.py | 13 ++-
trunk/src/mse/COVARIABLE/covariable.py | 91 +++++++++++++------
trunk/src/mse/MODELS/models.py | 61 +++++--------
trunk/src/mse/dynamicalSystem.py | 10 ++
trunk/tests/test_paramOpt.py | 2 +
6 files changed, 112 insertions(+), 76 deletions(-)
diff --git a/trunk/mseMatTools/Frontend/test_MseMatTools.py b/trunk/mseMatTools/Frontend/test_MseMatTools.py
index 34fe6c2..3efcf26 100644
--- a/trunk/mseMatTools/Frontend/test_MseMatTools.py
+++ b/trunk/mseMatTools/Frontend/test_MseMatTools.py
@@ -95,10 +95,13 @@ print(outArray)
outWait=np.array([1.0,2.0,3.0,4.0], dtype=np.double)
if (abs(LA.norm(outArray-outWait)) > 1e-10):
exit(2)
+outWait=np.array([-1.0,-2.0,-3.0,4.0], dtype=np.double)
+outArray=np.array([0.0,.0,.0,4.0], dtype=np.double)
+
outBuf=np.array([0.0,0.0,0.0,4.0], dtype=np.double)
MT.addInPlaceVec("testLoad",outArray,outBuf,-1)
-if (abs(LA.norm(outArray)) > 1e-10):
- exit(3)
-MT.addInPlaceVec("testLoad",outArray,outBuf,1)
if (abs(LA.norm(outArray-outWait)) > 1e-10):
- exit(4)
+ exit(3)
+
+def test_mat():
+ return 0
diff --git a/trunk/mseMatTools/Frontend/test_MseMatToolsRaster.py b/trunk/mseMatTools/Frontend/test_MseMatToolsRaster.py
index 5753546..a3f0933 100644
--- a/trunk/mseMatTools/Frontend/test_MseMatToolsRaster.py
+++ b/trunk/mseMatTools/Frontend/test_MseMatToolsRaster.py
@@ -1,5 +1,6 @@
import pymseMatTools.mseMatTools as MT
import numpy as np
+import os
NCOL=20
NROW=20
@@ -10,14 +11,16 @@ for i in range(20): Ycoords[i*NCOL:(i+1)*NCOL]=i
dofs=np.zeros(NCOL*NROW,np.double)
returnV=np.array([0],dtype=np.int32)
-MT.fillFromRaster("../RASTER/testRaster",dofs,Xcoords,Ycoords,returnV)
+absPath=os.getcwd().split("trunk")[0]+"/trunk/mseMatTools/RASTER/"
+MT.fillFromRaster(absPath+"/testRaster",dofs,Xcoords,Ycoords,returnV)
for i in range(20): print(dofs[i*NCOL:(i+1)*NCOL])
res=0
-for i in range(10):
+for i in range(10,20):
if (dofs[i*NCOL] != 1 or dofs[i*NCOL+5] != 1):
res=1
-for i in [NCOL-2,NCOL-1]:
+for i in [0,1]:
if (dofs[i*NCOL+NROW-1] != 1 or dofs[i*NCOL+NROW-4] != 1 or dofs[i*NCOL+NROW-5] != 0):
- res=1
+ res=2
print("test value="+str(res))
-exit(res)
+def test_raster():
+ return res
diff --git a/trunk/src/mse/COVARIABLE/covariable.py b/trunk/src/mse/COVARIABLE/covariable.py
index 02d5a01..6d27b44 100644
--- a/trunk/src/mse/COVARIABLE/covariable.py
+++ b/trunk/src/mse/COVARIABLE/covariable.py
@@ -52,32 +52,18 @@ class builderCovariableFromRaster(builderCovariable):
#
#
class covariable:
- def __init__(self,aDS,aName) -> None:
+ def __init__(self,aName,aDS) -> None:
self.dS=aDS
- self.dS.covariables.append(self)
self.name=aName
self.femCov=None
self.gradCov=None
- def save(self):
- MT.saveVec(self.dS.study.covDir+self.name,self.femCov.x.array)
- def load(self):
- MT.loadVec(self.dS.study.covDir+self.name,self.femCov.x.array)
- def exportVtk(self,fileName=None):
- f=None
- if (fileName):
- f=VTKFile(self.dS.dolfinMesh.comm, fileName+".pvd","w")
- else:
- f=VTKFile(self.dS.dolfinMesh.comm, self.dS.study.exportDir+self.name+".pvd","w")
- f.write_function(self.femCov)
- f.close()
-
class covariableFromExpression(covariable):
- def __init__(self,aDS,aName) -> None:
+ def __init__(self,aName,aDS) -> None:
super().__init__(aDS,aName)
self.femCov=Function(self.dS.vSpace)
self.femCov.x.array.fill(0.0)
self.gradCov=None
- def fillFromExpression(self,aExpression):
+ def setValue(self,aExpression):
def fct_express(x):
return eval(aExpression)
self.femCov.interpolate(fct_express)
@@ -89,26 +75,77 @@ class covariableFromExpression(covariable):
class covariableCst(covariable):
def __init__(self,aDS,aName) -> None:
super().__init__(aDS,aName)
- self.femCov=0
+ self.femCov=Constant(self.dS.dolfinMesh,ScalarType(0.0))
self.gradCov=0
def setValue(self,value):
- self.femCov=value
+ self.femCov.value(ScalarType(value))
def save(self):
pass
def load(self):
pass
def exportVtk(self,fileName=None):
pass
-class covariableFromRaster(covariable):
+class covariableFromCovFile(covariable):
def __init__(self,aDS,aName) -> None:
super().__init__(aDS,aName)
self.femCov=Function(self.dS.vSpace)
self.femCov.x.array.fill(0.0)
self.gradCov=None
- def fillFromRaster(self,rasterFileInfo):
- returnV=np.array([0],dtype=np.int32)
- XX= self.dS.vSpace.tabulate_dof_coordinates()[:,0]
- YY= self.dS.vSpace.tabulate_dof_coordinates()[:,1]
- MT.fillFromRaster(rasterFileInfo,self.femCov.x.array,XX,YY,returnV)
- if (returnV[0]):
- MSEPrint(" failed to import raster. err="+str(returnV[0]))
+ def exportVtk(self,fileName=None):
+ f=None
+ if (fileName):
+ f=VTKFile(self.dS.dolfinMesh.comm, fileName+".pvd","w")
+ else:
+ f=VTKFile(self.dS.dolfinMesh.comm, self.dS.study.exportDir+self.name+".pvd","w")
+ f.write_function(self.femCov)
+ f.close()
+ def setValue(self):
+ MT.loadVec(self.study.covDir+self.name,self.femCov.x.array)
+ #Pour les covariable , variable en temps, on creer un fichier timeCovName.txt contenant deux colonnes (temps, fichier covariable):
+ # covFileName0 0.0
+ # covFileName1 1.0
+ # covFileName2 1.4
+ # covFileName3 2.0
+ # par exemple pour les donnée de températures.
+ # A l'instar de l'explorTraj, la mise a jours de la covariable pendant la simulation, constera à chercher l'interval de temp correspondant au pas de simu et a faire l'interpolation, (ou le pas précédent)
+ def setValue(self,aTime):
+ prevFile=""
+ curFile=""
+
+ f_Time = open(self.study.covDir+self.name+"Time.txt")
+ if (not f_time):
+ MSEPrint(" Erreur, no time file for covariable "+self.study.covDir+self.name+"Time.txt")
+ return
+ line = f_Time.readline()
+ #print("line = ",line)
+ curT=float(line.split()[0])
+ prevT=curT
+ curFile=int(line.split()[1])
+ while line := f_Time.readline():
+ prevT=curT
+ prevFile=nextFile
+ curT=float(line.split()[0])
+ curFile=int(line.split()[1])
+ if ((curT-aTime>-1e-16) == (curT-prevT>-1e-16)):
+ break
+ f_Time.close()
+ aCoef=0
+ if (abs(aTime-prevT) < 1e-16):
+ MT.loadVec(prevFile,self.femCov.x.array)
+ elif(abs(curT-aTime)<1e-16):
+ MT.loadVec(curFile,self.femCov.x.array)
+ else:
+ aCoef=(aTime-prevT)/(curT-prevT)
+ #print(str(aCoef))
+ if (aCoef <0 or aCoef > 1):
+ if (aCoef>1):
+ aCoef=1
+ else:
+ aCoef=0
+ for aDs in self.study.system.dynamicalSystems:
+ MT.loadVec(prevFile,self.dS.bufState.comps[0].x.array)
+ #use self.dS.bufState.comps[0] as buffer
+ aDs.trajState.load(self.trajDirName+"ds"+str(aDs.id),prevStep)
+ self.dS.bufState.comps[0].x.array[]=(1-aCoef)*self.dS.bufState.comps[0].x.array[]
+ MT.loadVec(curFile,self.femCov.x.array)
+ self.femCov.x.array[]=aCoef*self.femCov.x.array[]+self.dS.bufState.comps[0].x.array[]
diff --git a/trunk/src/mse/MODELS/models.py b/trunk/src/mse/MODELS/models.py
index f43ada2..87bd215 100644
--- a/trunk/src/mse/MODELS/models.py
+++ b/trunk/src/mse/MODELS/models.py
@@ -25,6 +25,23 @@ def strForDolf(strFile, unknownName="abcdefgh"):
"""
return(str(sympify(strFile).subs(list(zip(symbols(' '.join(unknownName)), symbols(' '.join(['$'+k for k in unknownName])))))))
+from re import sub
+
+# Fonction add sufix and prefix to upcase words ie covariable
+#suffix = "]"
+#prefix = "self.dS.dicCov["
+#nouveau_texte = ajouter_suffixe_majuscule(texte, prefixe, suffixe)
+#print(nouveau_texte)
+def formatCovariables(chaine, prefixe, suffixe):
+ # L'expression reguliere cherche des mots entierement en majuscules
+ regex = r'\b[A-Z]+\b'
+ # Fonction de remplacement qui ajoute le suffixe
+ def ajouter_suffixe(match):
+ return prefixe + match.group(0) + suffixe
+ # Utilise re.sub() pour remplacer les mots en majuscules par eux-memes avec le suffixe ajoute
+ return sub(regex, ajouter_suffixe, chaine)
+
+
class model:
#chose your Debug mode
#myPrint=MSEPrint
@@ -116,43 +133,6 @@ class model:
def getFtm1(self):
return self.Ftm1
-separateurs = [' ', ',', ';', '.',]
-
-def is_separator(char):
- return char in " ,;.:/*-+={([\^)]}"
-
-def replaceWordbyWord( inStr, oldW, newW):
- lDebug=0
- if (lDebug): model.myPrint("instr= "+inStr)
- prevIndex=0
- curIndex=inStr.find(oldW)
- outStr=""
- lenInStr=len(inStr)
- lenOldW=len(oldW)
- if (lDebug): model.myPrint("lenInStr="+str(lenInStr))
- if (lDebug): model.myPrint("lenOldW="+str(lenOldW))
- while (curIndex>=0):
- outStr=outStr+inStr[prevIndex:curIndex]
- if (lDebug): model.myPrint("curIndex="+str(curIndex))
- if (lDebug): model.myPrint("curIndex+lenOldW-1="+str(curIndex+lenOldW))
- if ((curIndex==0 or is_separator(inStr[curIndex-1])) and (curIndex+lenOldW == lenInStr or is_separator(inStr[curIndex+lenOldW]))):
- outStr=outStr+newW
- else:
- outStr=outStr+oldW
- prevIndex=curIndex+lenOldW
- if (lDebug): model.myPrint("prevIndex="+str(prevIndex))
- curIndex=inStr[prevIndex:].find(oldW)
- curIndex=len(inStr)
- outStr=outStr+inStr[prevIndex:curIndex]
- if (lDebug): model.myPrint("outStr= "+outStr)
- return outStr
-
-def updadeVarfWithCovariables(inStr,covs):
- i=0
- for cov in covs:
- inStr=replaceWordbyWord(inStr, cov.name, "self.dS.covariables["+str(i)+"].femCov")
- i+=1
- return inStr
class txtModel(model):
"""short presentation of class
@@ -216,7 +196,7 @@ class txtModel(model):
#for the comp number numComp \in {0,..,self.nbComps-1}:
#self.varfs[numComp*(self.nbComps+1)] contains the varf F_numComp for the rhs
#self.varfs[numComp*(self.nbComps+1)+(j+1)] contains the derivative d_jF_numComp, j \in {0,..,self.nbComps-1}
- strkeys=r"[\s\+\-\*/\=\(\)\^]+"
+ #strkeys=r"[\s\+\-\*/\=\(\)\^]+"
for numComps in range(self.nbComps):
##READ SOURCE TERM
strline=f.readline().split('=')[1].strip()
@@ -227,7 +207,7 @@ class txtModel(model):
strline=strline.replace('p'+str(count+1),'self.dS.paramD['+str(count)+']')
model.myPrint(" read txt model F"+txtModel.unknownName[numComps]+"="+strline)
#### replace COVARIABLES
- strline=updadeVarfWithCovariables(strline,self.dS.covariables)
+ strline=formatCovariables(strline, "self.dS.dicCov[\"", "\"]")
model.myPrint(strline)
#### BUILD VARF USING FK
code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+"*self.dS.v*self.dS.dx)"
@@ -244,9 +224,10 @@ class txtModel(model):
##replace param
for count,p in enumerate(self.dS.paramD):
strline=strline.replace('p'+str(count+1),'self.dS.paramD['+str(count)+']')
- strline=updadeVarfWithCovariables(strline,self.dS.covariables)
+ strline=formatCovariables(strline, "self.dS.dicCov[\"", "\"]")
model.myPrint("txtModel: read txt model d"+txtModel.unknownName[k]+"F"+txtModel.unknownName[numComps]+"="+strline)
code="self.varfs["+str(numComps)+"*(self.nbComps+1)+"+str(k)+"+1]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
+ model.myPrint("code for varf "+code)
exec(compile(code, 'sumstring', 'exec'))
aDs.addSourcesTerm(self)
def updateBlockMat(self,compi,compj):
diff --git a/trunk/src/mse/dynamicalSystem.py b/trunk/src/mse/dynamicalSystem.py
index c2b897e..361d456 100644
--- a/trunk/src/mse/dynamicalSystem.py
+++ b/trunk/src/mse/dynamicalSystem.py
@@ -79,8 +79,18 @@ class dynamicalSystem:
self.massMat=matrix(assemble_matrix(form(aux)))
self.covariables=[]
self.study.system._addDS(self)
+ #diconary of covariables
+ self.dicHeterogenTimeVarCov = dict()
+ self.dicHeterogenCov = dict()
+ self.dicHomogenCov = dict()
#trajDirectory
#self.trajState = explorTrajectory(self.study, self.trajState)
+ def addHomogenCov(self,aName,aCov):
+ self.dicHomogenCov[aName]=aCov
+ def addHeterogenTimeVarCov(self,aName,aCov):
+ self.dicHeterogenTimeVarCov[aName]=aCov
+ def addHeterogenCov(self,aName,aCov):
+ self.dicHeterogenCov[aName]=aCov
def setParamDolf(self):
"""updating the param with the study param
"""
diff --git a/trunk/tests/test_paramOpt.py b/trunk/tests/test_paramOpt.py
index 14451fb..a48c6a2 100644
--- a/trunk/tests/test_paramOpt.py
+++ b/trunk/tests/test_paramOpt.py
@@ -100,6 +100,8 @@ def aCallableOpt(intermediate_result: OptimizeResult):
def test_paramOpt(defaultParam):
""" call function to find param opt
"""
+ #to long
+ #return 0
global globalCount
# SIMU #
#D = 5.5
--
GitLab
From 9ce88361ab5104e42792157c516bedbaabafb5b6 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Thu, 30 May 2024 11:15:49 +0200
Subject: [PATCH 08/23] move test_paramOpt to DEMOS/study8.py
---
.../test_paramOpt.py => DEMOS/study8.py} | 26 +++++++------------
1 file changed, 10 insertions(+), 16 deletions(-)
rename trunk/{tests/test_paramOpt.py => DEMOS/study8.py} (86%)
diff --git a/trunk/tests/test_paramOpt.py b/trunk/DEMOS/study8.py
similarity index 86%
rename from trunk/tests/test_paramOpt.py
rename to trunk/DEMOS/study8.py
index a48c6a2..938083c 100644
--- a/trunk/tests/test_paramOpt.py
+++ b/trunk/DEMOS/study8.py
@@ -15,7 +15,6 @@ from os.path import isdir, exists
import numpy as np
import matplotlib.pyplot as plt
import shutil
-import pytest
from scipy.optimize import minimize, OptimizeResult
@@ -33,7 +32,7 @@ testNLTS=1
#une etude
-studName="TEST_PARAMOPT"
+studName="STUDY8"
aStudy=study(studName, nbParam=2)
aStudy.setGeometry("SQUARE",0.4)
@@ -93,27 +92,20 @@ globalCount=0
def aCallableOpt(intermediate_result: OptimizeResult):
"""a callable use in intermediate odo a kdf ozf zkf sofnqqk oajn
"""
- print("val of param in init: ",intermediate_result.x)
- paramPath.append(intermediate_result.x)
+ print("In iteration: \n",intermediate_result)
+ paramPath.append(intermediate_result.x.copy())
-@pytest.mark.parametrize("defaultParam",[([0.9,1]),([1.2,1.]),([1,0.54]),([1,1.3]),([1.4,1.3]),([0.64,0.73])])
-def test_paramOpt(defaultParam):
+def findParamOpt(defaultParam):
""" call function to find param opt
"""
- #to long
- #return 0
global globalCount
# SIMU #
- #D = 5.5
- #aDiffModel.setDiffCoef(D)
print("********\nParam a_sim =",defaultParam)
aStudy.setParam(defaultParam)
aStudy.simulator.doSimu()
#we call function to find paramOpt
- optParam = minimize(quaDiff.sysAdj.valJAndGradJ, defaultParam, method="L-BFGS-B", bounds=[(0.1,3)], jac=True, tol=None, callback=aCallableOpt)
+ optParam = minimize(quaDiff.sysAdj.valJAndGradJ, defaultParam, method="L-BFGS-B", bounds=[(0.1,3)], jac=True, tol=None, callback=aCallableOpt, options={"maxiter":7})
print (optParam)
- print(optParam.nit)
- print (optParam.x)
tmpParamPath = paramPath[-1-optParam.nit:]
tmpParamPath.insert(0,np.array(defaultParam))
tmpParamPath = np.array(tmpParamPath)
@@ -128,6 +120,8 @@ def test_paramOpt(defaultParam):
plt.axis((0.5,1.5,0.5,1.5))
globalCount+=1
plt.title("parameter for test "+str(globalCount))
- plt.savefig("pathParameter"+str(globalCount)+".png")
- assert(abs(optParam.x[0]-1.)<2e-4)
- assert(abs(optParam.x[1]-1.)<2e-4)
+ plt.savefig(studName+"/pathParameter"+str(globalCount)+".png")
+### list of param to test
+paramInit = [([0.9,1]),([1.2,1.]),([1,0.54]),([1,1.3]),([1.4,1.3]),([0.64,0.73])]
+### call of function to find paramOPt
+findParamOpt(paramInit[0])
--
GitLab
From a48ebe208701e6916e61e851deb5d24c677d90ee Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Thu, 30 May 2024 11:18:25 +0200
Subject: [PATCH 09/23] simplify test: finite difference forward instead of FD
central, less param init => less test
---
trunk/tests/test_backwardModel.py | 23 ++++++++++++-----------
1 file changed, 12 insertions(+), 11 deletions(-)
diff --git a/trunk/tests/test_backwardModel.py b/trunk/tests/test_backwardModel.py
index 43f4755..9e6af36 100644
--- a/trunk/tests/test_backwardModel.py
+++ b/trunk/tests/test_backwardModel.py
@@ -90,7 +90,7 @@ quaDiff=quadraticDiff(ds2d1)
# We test with different epsilon and parameters.
print("Param a_ref =",aStudy.getParam)
-@pytest.mark.parametrize("defaultParam",[([0.5,0.1]),([1.5,1.]),([0.9,0.54]),([1.2,0.6])])
+@pytest.mark.parametrize("defaultParam",[([0.5,0.1]),([1.5,1.])])
def test_GradWithFiniteDifference(defaultParam):
"""
test the approximation of GradV mse with finite Difference.
@@ -105,29 +105,30 @@ def test_GradWithFiniteDifference(defaultParam):
aStudy.setParam(defaultParam)
aStudy.simulator.doSimu()
#
+ VIn = quaDiff.computeV()
quaDiff.computeGradV()
GJ = quaDiff.gradJ
print("gradV = :",GJ)
param= [None]*aStudy.nbParam
gradDif = [None]*aStudy.nbParam
- tabEpsilon = [1e-2, 1e-3]
+ tabEpsilon = [5e-3, 5e-4]
for eps in tabEpsilon:
print("eps for diff finie = ",eps)
invEps = 1/eps
for i in range(aStudy.nbParam):
#define param
- param= [defaultParam[k] + (k==i)*eps*0.5 for k in range(aStudy.nbParam)]
+ param= [defaultParam[k] + (k==i)*eps for k in range(aStudy.nbParam)]
aStudy.setParam(param)
#Simu
aStudy.simulator.doSimu()
VForeward = quaDiff.computeV()
- #define param
- param= [defaultParam[k] - (k==i)*eps*0.5 for k in range(aStudy.nbParam)]
- aStudy.setParam(param)
- #Simu
- aStudy.simulator.doSimu()
- #val of V
- gradDif[i] = (VForeward - quaDiff.computeV())*invEps
+ ##define param
+ #param= [defaultParam[k] - (k==i)*eps*0.5 for k in range(aStudy.nbParam)]
+ #aStudy.setParam(param)
+ ##Simu
+ #aStudy.simulator.doSimu()
+ ##val of V
+ gradDif[i] = (VForeward - VIn)*invEps
print("gradV of finite element = ",gradDif)
print("gradV of MSE = ", GJ)
print("dir grad FE: ",gradDif[1]*1./gradDif[0])
@@ -144,7 +145,7 @@ def test_GradWithFiniteDifference(defaultParam):
print("Erreur direction: "+str(errDir)+"%")
print("Erreur norm2: "+str(errNorm)+"%")
assert(errNorm < 1.5*1e-1)
- assert(errDir < 3.5*1e-1)
+ assert(errDir < 4.*1e-1)
#for i in range(aStudy.nbParam):
# tmpscal = gradDif[i] if gradDif[i] !=0 else 1
# assert(abs((gradDif[i] - GJ[i])/tmpscal)< 5*1e-2)
--
GitLab
From f9a4d62a046a323fc3a6443d95e6c76c84893810 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Thu, 30 May 2024 11:44:49 +0200
Subject: [PATCH 10/23] Patch init of obsProcess + in test, remove call of
computeGradV, we only have unit tests
---
trunk/src/mse/obsProcess.py | 4 ++--
trunk/tests/test_obsProcess.py | 24 +++++++-----------------
2 files changed, 9 insertions(+), 19 deletions(-)
diff --git a/trunk/src/mse/obsProcess.py b/trunk/src/mse/obsProcess.py
index 5e62f76..39481f9 100644
--- a/trunk/src/mse/obsProcess.py
+++ b/trunk/src/mse/obsProcess.py
@@ -21,7 +21,7 @@ class obsProcess:
This class contain the definition of main function for an obsprocess.
"""
- def __init__(self, aDS, aDirURef="UREF/")->None:
+ def __init__(self, aDS, aDirURef=None)->None:
#Dynamical system
self.dS=aDS
self.dS.setObsProcess(self)
@@ -34,7 +34,7 @@ class obsProcess:
# array containing grad of V (J), length nbParam
self.gradJ=np.zeros(self.dS.study.nbParam)
#directory containing file with value of uRef
- if (aDirURef==""):
+ if (aDirURef==None):
self.dirURef=self.dS.study.absoluteStudyDir+"/UREF/"
else:
self.dirURef=aDirURef
diff --git a/trunk/tests/test_obsProcess.py b/trunk/tests/test_obsProcess.py
index b0428c3..3929b1f 100644
--- a/trunk/tests/test_obsProcess.py
+++ b/trunk/tests/test_obsProcess.py
@@ -93,22 +93,12 @@ rename(path+"tmptime.txt", path+"time.txt")
aStudy.setParam([2.,4.])
aStudy.simulator.doSimu()
-def test_ClassObsprocess():
- """test creation of class obsProcess
+def test_ClassObsProcess_Ds():
+ """test of class obsProcess and dynamicalSystem
"""
- assert True
-
-
-def test_FnctObsComputeGradV():
- """test function computeGradV from class obsProcess
-
- we test first if file are create and second if gradV is correct
+ assert (obsP.dS == ds2d1)
+ assert(ds2d1.obsProcess == obsP)
+def test_ClassObsProcess_Dir():
+ """test if directory is correct in class obsProcess
"""
- obsP.computeGradV()
- #test file
- for sourceTerm in obsP.dS.sourcesTerm:
- assert (os.path.exists(sourceTerm.txtFile+"Dp.txt"))
- for i in range(len(obsP.gradJ)):
- print("gradJ["+str(i)+"] = ",obsP.gradJ[i])
- #test value
- assert True
+ assert (obsP.dirURef == aStudy.absoluteStudyDir+"/UREF/")
--
GitLab
From f7d33db37e0c58776abfa3d0a3971ba28f4d2fcb Mon Sep 17 00:00:00 2001
From: Olivier Bonnefon <olivier.bonnefon@inra.fr>
Date: Fri, 31 May 2024 10:09:40 +0200
Subject: [PATCH 11/23] fix bug in obsProcess, cancel verbose mode in
mseMatTools, update README.md
---
trunk/README.md | 1 -
trunk/mseMatTools/src/mseMatTools.c | 2 +-
trunk/src/mse/quadraticDiff.py | 2 +-
3 files changed, 2 insertions(+), 3 deletions(-)
diff --git a/trunk/README.md b/trunk/README.md
index 0c3eac9..2c33d91 100644
--- a/trunk/README.md
+++ b/trunk/README.md
@@ -37,7 +37,6 @@ echo "remove previuos builds"
rm -rf mseMatTools/build/* mseMatTools/Frontend/build/* mseMapTools/build/* mseMapTools/Frontend/build/*
echo "1) install mseMatTools :"
-guix time-machine -C ../guix/myCurrentChannels.scm -- shell -C --share=/tmp --preserve='^DISPLAY$' -m ../guix/manifestCompilmseMatTools.scm -m ../guix/manifestCompilpymseMatTools.scm -m ../guix/manifestToCompilMSE.scm -m ../guix/manifestDolfinxv06.scm
cd $TRUNK_MSE
mkdir mseMatTools/build
diff --git a/trunk/mseMatTools/src/mseMatTools.c b/trunk/mseMatTools/src/mseMatTools.c
index 8276756..edd51c6 100644
--- a/trunk/mseMatTools/src/mseMatTools.c
+++ b/trunk/mseMatTools/src/mseMatTools.c
@@ -2,7 +2,7 @@
#include <petscmat.h>
#include "stdio.h"
#include <stdlib.h>
-#define mseMatToolsS_VERBOSE
+//#define mseMatToolsS_VERBOSE
void getPetscId(Mat aMat){
printf("Begin getPetscId2\n");
diff --git a/trunk/src/mse/quadraticDiff.py b/trunk/src/mse/quadraticDiff.py
index 8a5f030..fe1fce3 100644
--- a/trunk/src/mse/quadraticDiff.py
+++ b/trunk/src/mse/quadraticDiff.py
@@ -25,7 +25,7 @@ class quadraticDiff(obsProcess):
:param str aDirURef: an absolute pasth to the directory where observation files are. By default in dir study.absoulteSutyDir+'UREF/'.
"""
myPrint=MSEPrint_0
- def __init__(self,aDS, aDirURef="") -> None:
+ def __init__(self,aDS, aDirURef=None) -> None:
super().__init__(aDS, aDirURef)
#explor_uRef type of explorTrajectory
self.explor_uRef = explorTrajectory(self.dS.study, self.dirURef)
--
GitLab
From 5c06eb81c1930886c8144bde32698c7cd3540af5 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Fri, 31 May 2024 15:26:30 +0200
Subject: [PATCH 12/23] update init of txtDiffusionModel, not finish
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 165 ++++++++++----------
1 file changed, 86 insertions(+), 79 deletions(-)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 60ec276..4daaa63 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -9,7 +9,6 @@ from mse.linAlgTools import matrix,vector
class diffusionModel(model):
def __init__(self,aDs,aDiffCoef=1) -> None:
super().__init__(aDs)
- ### NOTE Test: we define diffusionModel with param
aDs.addSecondOrderTerm(self)
self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
@@ -43,95 +42,103 @@ class txtDiffusionModel(model):
"""
def __init__(self,aDs, aTxtFile) -> None:
super().__init__(aDs)
+ aDs.addSecondOrderTerm(self)
+ self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
+ self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
+ self.block=None
+ self.blockAdjoint=None
+ self.bilinearForm=None
+ self._computeBlock()
+ self.updateLinearisedMat()
+ def setDiffCoef(self,aTxtFile):
+ self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
+ self._computeBlock()
+ self.updateLinearisedMat()
+ def _computeBlock(self):
+ aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
+ self.bilinearForm=form(aux)
+ self.block=matrix(assemble_matrix(self.bilinearForm))
# we check if file exist
- if (exists(aTxtFile+".txt")):
- self.txtFile=aTxtFile
- else:
- self.txtFile=computerEnv.modelsDir+aTxtFile #aTxtFile form is "affineGrowth"
- model.myPrint("Looking for model "+self.txtFile)
- try:
- f = open(self.txtFile+".txt", "r")
- except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+".txt")
- return
- # we check if the derivate of parameter file exist
- try:
- fd = open(self.txtFile+"DVP.txt", "r")
- except FileNotFoundError:
- # if didn't exist, we create it
- self.dvFileExist=False
- print("Création fichier "+self.txtFile+"DVP.tx")
- fd = open(self.txtFile+"DVP.txt", "w") #open file in mod write
- for numComps in range(self.nbComps):
- strF=f.readline()
- strF=strF.split('=')[1].strip()#.replace("^","**")
- for p in range(self.aDs.study.nbParam):
- grad1s=diff(strF,'p'+str(p+1))
- grad1=str(grad1s).replace("**","^")
- #write in fic
- fd.write("dp"+str(p+1)+"f"+txtModel.unknownName[numComps]+"="+grad1+"\n")
- f.close()
- fd.close()
- f = open(self.txtFile+".txt", "r")
+ if (exists(aTxtFile+".txt")):
+ self.txtFile=aTxtFile
+ else:
+ self.txtFile=computerEnv.modelsDir+aTxtFile #aTxtFile form is "affineGrowth"
+ model.myPrint("Looking for model "+self.txtFile)
+ try:
+ f = open(self.txtFile+".txt", "r")
+ except IOError:
+ model.myPrint("IOError in txtModel when reading "+self.txtFile+".txt")
+ return
+ # we check if the derivate of parameter file exist
try:
fd = open(self.txtFile+"DVP.txt", "r")
+ except FileNotFoundError:
+ # if didn't exist, we create it
+ self.dvFileExist=False
+ print("Création fichier "+self.txtFile+"DVP.tx")
+ fd = open(self.txtFile+"DVP.txt", "w") #open file in mod write
+ for numComps in range(self.nbComps):
+ strF=f.readline()
+ strF=strF.split('=')[1].strip()#.replace("^","**")
+ for p in range(self.aDs.study.nbParam):
+ grad1s=diff(strF,'p'+str(p+1))
+ grad1=str(grad1s).replace("**","^")
+ #write in fic
+ fd.write("dp"+str(p+1)+"f"+txtModel.unknownName[numComps]+"="+grad1+"\n")
+ f.close()
+ fd.close()
+ f = open(self.txtFile+".txt", "r")
+ try:
+ fd = open(self.txtFile+"DVP.txt", "r")
+ except IOError:
+ model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
+ return
except IOError:
model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
return
- except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
- return
- ### find idea for derivate in obprocess.computeGradV
- ### TODO: read file + create varfs like txtmodel
- ###Â we have \Delta H(u),
- ### H_comps, Dp_{ai}H_comps
- self.varfs=[None]*(self.nbComps*(self.dS.study.nbParam+1))
- self.varfsTm1=[None]*self.nbComps
- #for the comp number numComp \in {0,..,self.nbComps-1}:
- #self.varfs[numComp*(self.nbComps+1)] contains the varf F_numComp for the rhs
- #self.varfs[numComp*(self.nbComps+1)+(j+1)] contains the derivative d_jF_numComp, j \in {0,..,self.nbComps-1}
- strkeys=r"[\s\+\-\*/\=\(\)\^]+"
- for numComps in range(self.nbComps):
- ##READ SECOND ORDER TERM
- strline=f.readline().split('=')[1].strip()
- strline=strForDolf(strline, self.unknownName)
- strline=strline.replace('^','**')
- ##replace param
- for count_p in range(self.dS.paramD):
- strline=strline.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
- model.myPrint(" read txt model F"+txtModel.unknownName[numComps]+"="+strline)
- #### replace COVARIABLES
- strline=updadeVarfWithCovariables(strline,self.dS.covariables)
- model.myPrint(strline)
- #### BUILD VARF USING FK
- code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+"*self.dS.v*self.dS.dx)"
- exec(compile(code, 'sumstring', 'exec'))
- #### BUILD VARF USING UTM1
- code="self.varfsTm1["+str(numComps)+"]="+self._subFXForm(strline,"utm1")+"*self.dS.v*self.dS.dx)"
- model.myPrint(code)
- exec(compile(code, 'sumstring', 'exec'))
- ## READ DERIVATIVES PARAM
- for countP in range(self.dS.paramD):
- strline=fd.readline().split('=')[1].strip()
- strline=strForDolf(strline,self.unknownName)
- #strline=strline.replace('^','**')
+ ### find idea for derivate in obprocess.computeGradV
+ ### TODO: read file + create varfs like txtmodel
+ ###Â we have \Delta H(u),
+ ### H_comps, Dp_{ai}H_comps
+ self.varfs=[None]*(self.nbComps*(self.dS.study.nbParam+1))
+ self.varfsTm1=[None]*self.nbComps
+ #for the comp number numComp \in {0,..,self.nbComps-1}:
+ #self.varfs[numComp*(self.nbComps+1)] contains the varf F_numComp for the rhs
+ #self.varfs[numComp*(self.nbComps+1)+(j+1)] contains the derivative d_jF_numComp, j \in {0,..,self.nbComps-1}
+ strkeys=r"[\s\+\-\*/\=\(\)\^]+"
+ for numComps in range(self.nbComps):
+ ##READ SECOND ORDER TERM
+ strline=f.readline().split('=')[1].strip()
+ strline=strForDolf(strline, self.unknownName)
+ strline=strline.replace('^','**')
##replace param
- for tmpCountP in enumerate(self.dS.paramD):
- strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountp)+']')
+ for count_p in range(self.dS.paramD):
+ strline=strline.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
+ model.myPrint(" read txt model F"+txtModel.unknownName[numComps]+"="+strline)
+ #### replace COVARIABLES
strline=updadeVarfWithCovariables(strline,self.dS.covariables)
- model.myPrint("txtModel: read txt model d"+txtModel.unknownName[k]+"F"+txtModel.unknownName[numComps]+"="+strline)
- code="self.varfs["+str(numComps)+"*(self.dS.study.nbParam+1)+"+str(k)+"+1]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
+ model.myPrint(strline)
+ #### BUILD VARF USING FK
+ code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+"*self.dS.v*self.dS.dx)"
exec(compile(code, 'sumstring', 'exec'))
- def _computeBlock(self):
- aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
- self.bilinearForm=form(aux)
- self.block=matrix(assemble_matrix(self.bilinearForm))
+ #### BUILD VARF USING UTM1
+ code="self.varfsTm1["+str(numComps)+"]="+self._subFXForm(strline,"utm1")+"*self.dS.v*self.dS.dx)"
+ model.myPrint(code)
+ exec(compile(code, 'sumstring', 'exec'))
+ ## READ DERIVATIVES PARAM
+ for countP in range(self.dS.paramD):
+ strline=fd.readline().split('=')[1].strip()
+ strline=strForDolf(strline,self.unknownName)
+ #strline=strline.replace('^','**')
+ ##replace param
+ for tmpCountP in enumerate(self.dS.paramD):
+ strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountp)+']')
+ strline=updadeVarfWithCovariables(strline,self.dS.covariables)
+ model.myPrint("txtModel: read txt model d"+txtModel.unknownName[k]+"F"+txtModel.unknownName[numComps]+"="+strline)
+ code="self.varfs["+str(numComps)+"*(self.dS.study.nbParam+1)+"+str(k)+"+1]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
+ exec(compile(code, 'sumstring', 'exec'))
def updateBlockMat(self,compi,compj):
return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)+compj+1]))
- def _computeBlockAdjoint(self):
- aux=dot(grad(self.dS.u),grad(self.diffCoef*self.dS.v))*self.dS.dx
- self.bilinearForm=form(aux)
- self.blockAdjoint=matrix(assemble_matrix(self.bilinearForm))
def updateBlockMatAdjoint(self,compi,compj):
if (compi==compj):
return self.blockAdjoint
--
GitLab
From 03641614351c3bb163ac0fc8bde1c087b7390006 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Fri, 31 May 2024 15:34:21 +0200
Subject: [PATCH 13/23] Deletes modif in diffuisionModel that added a
parameter. NOTE: this modif will be moved to another class
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 4 ++++
trunk/src/mse/obsProcess.py | 2 +-
2 files changed, 5 insertions(+), 1 deletion(-)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 4daaa63..7cbb4f9 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -9,8 +9,12 @@ from mse.linAlgTools import matrix,vector
class diffusionModel(model):
def __init__(self,aDs,aDiffCoef=1) -> None:
super().__init__(aDs)
+<<<<<<< HEAD
aDs.addSecondOrderTerm(self)
self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
+=======
+ self.coefName=aDiffCoef
+>>>>>>> 57c81bb (Deletes modif in diffuisionModel that added a parameter. NOTE: this modif will be moved to another class)
self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
self.block=None
self.blockAdjoint=None
diff --git a/trunk/src/mse/obsProcess.py b/trunk/src/mse/obsProcess.py
index 39481f9..5b3e8dd 100644
--- a/trunk/src/mse/obsProcess.py
+++ b/trunk/src/mse/obsProcess.py
@@ -199,7 +199,7 @@ class obsProcess:
#read u_a
code="self.daF.x.array[:]="+'+'.join(strdFp[i])
exec(compile(code,'sumstring', 'exec'))
- curF[i] = assemble_scalar(form(comp*self.daF*self.dS.dx)) + assemble_scalar(form(-inner(grad(self.dS.trajState.comps[count]),grad(comp))*self.dS.dx))*(self.dS.study.param[i]==self.dS.secondOrderTerm.coefName)
+ curF[i] = assemble_scalar(form(comp*self.daF*self.dS.dx))
### End V1
### V2
# NOTE: cette version demande de creer un vecteur. L'impementation est à l'abandon
--
GitLab
From 218d0832d8d27491208cb883e426c5ffe68319a6 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Thu, 30 May 2024 10:25:09 +0200
Subject: [PATCH 14/23] quick update, not finish
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 30 +++++++++++++++++++++
1 file changed, 30 insertions(+)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 7cbb4f9..a40920b 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -9,12 +9,21 @@ from mse.linAlgTools import matrix,vector
class diffusionModel(model):
def __init__(self,aDs,aDiffCoef=1) -> None:
super().__init__(aDs)
+<<<<<<< HEAD
<<<<<<< HEAD
aDs.addSecondOrderTerm(self)
self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
=======
self.coefName=aDiffCoef
>>>>>>> 57c81bb (Deletes modif in diffuisionModel that added a parameter. NOTE: this modif will be moved to another class)
+=======
+ self.coefName=aDiffCoef
+=======
+ ### NOTE Test: we define diffusionModel with param
+ aDs.addSecondOrderTerm(self)
+ self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
+>>>>>>> 3f748d9 (quick update, not finish)
+>>>>>>> 19312c4 (quick update, not finish)
self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
self.block=None
self.blockAdjoint=None
@@ -100,6 +109,7 @@ class txtDiffusionModel(model):
except IOError:
model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
return
+<<<<<<< HEAD
### find idea for derivate in obprocess.computeGradV
### TODO: read file + create varfs like txtmodel
###Â we have \Delta H(u),
@@ -143,6 +153,26 @@ class txtDiffusionModel(model):
exec(compile(code, 'sumstring', 'exec'))
def updateBlockMat(self,compi,compj):
return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)+compj+1]))
+=======
+ except IOError:
+ model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
+ return
+ ### find idea for derivate in obprocess.computeGradV
+ ### TODO: read file + create varfs like txtmodel
+ def truc():
+ pass
+ def _computeBlock(self):
+ aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
+ self.bilinearForm=form(aux)
+ self.block=matrix(assemble_matrix(self.bilinearForm))
+ def updateBlockMat(self,compi,compj):
+ if (compi==compj):
+ return self.block
+ def _computeBlockAdjoint(self):
+ aux=dot(grad(self.dS.u),grad(self.diffCoef*self.dS.v))*self.dS.dx
+ self.bilinearForm=form(aux)
+ self.blockAdjoint=matrix(assemble_matrix(self.bilinearForm))
+>>>>>>> 19312c4 (quick update, not finish)
def updateBlockMatAdjoint(self,compi,compj):
if (compi==compj):
return self.blockAdjoint
--
GitLab
From 80e49d180fb3a86df856d9784e2e175035fd495b Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Fri, 31 May 2024 15:26:30 +0200
Subject: [PATCH 15/23] update init of txtDiffusionModel, not finish
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 91 +++++++++++++++++++++
1 file changed, 91 insertions(+)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index a40920b..b816abc 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -10,6 +10,7 @@ class diffusionModel(model):
def __init__(self,aDs,aDiffCoef=1) -> None:
super().__init__(aDs)
<<<<<<< HEAD
+<<<<<<< HEAD
<<<<<<< HEAD
aDs.addSecondOrderTerm(self)
self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
@@ -17,9 +18,13 @@ class diffusionModel(model):
self.coefName=aDiffCoef
>>>>>>> 57c81bb (Deletes modif in diffuisionModel that added a parameter. NOTE: this modif will be moved to another class)
=======
+=======
+>>>>>>> ef8e663 (update init of txtDiffusionModel, not finish)
self.coefName=aDiffCoef
=======
### NOTE Test: we define diffusionModel with param
+=======
+>>>>>>> 5c06eb8 (update init of txtDiffusionModel, not finish)
aDs.addSecondOrderTerm(self)
self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
>>>>>>> 3f748d9 (quick update, not finish)
@@ -67,6 +72,7 @@ class txtDiffusionModel(model):
self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
self._computeBlock()
self.updateLinearisedMat()
+<<<<<<< HEAD
def _computeBlock(self):
aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
self.bilinearForm=form(aux)
@@ -161,11 +167,93 @@ class txtDiffusionModel(model):
### TODO: read file + create varfs like txtmodel
def truc():
pass
+=======
+>>>>>>> ef8e663 (update init of txtDiffusionModel, not finish)
def _computeBlock(self):
aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
self.bilinearForm=form(aux)
self.block=matrix(assemble_matrix(self.bilinearForm))
+ # we check if file exist
+ if (exists(aTxtFile+".txt")):
+ self.txtFile=aTxtFile
+ else:
+ self.txtFile=computerEnv.modelsDir+aTxtFile #aTxtFile form is "affineGrowth"
+ model.myPrint("Looking for model "+self.txtFile)
+ try:
+ f = open(self.txtFile+".txt", "r")
+ except IOError:
+ model.myPrint("IOError in txtModel when reading "+self.txtFile+".txt")
+ return
+ # we check if the derivate of parameter file exist
+ try:
+ fd = open(self.txtFile+"DVP.txt", "r")
+ except FileNotFoundError:
+ # if didn't exist, we create it
+ self.dvFileExist=False
+ print("Création fichier "+self.txtFile+"DVP.tx")
+ fd = open(self.txtFile+"DVP.txt", "w") #open file in mod write
+ for numComps in range(self.nbComps):
+ strF=f.readline()
+ strF=strF.split('=')[1].strip()#.replace("^","**")
+ for p in range(self.aDs.study.nbParam):
+ grad1s=diff(strF,'p'+str(p+1))
+ grad1=str(grad1s).replace("**","^")
+ #write in fic
+ fd.write("dp"+str(p+1)+"f"+txtModel.unknownName[numComps]+"="+grad1+"\n")
+ f.close()
+ fd.close()
+ f = open(self.txtFile+".txt", "r")
+ try:
+ fd = open(self.txtFile+"DVP.txt", "r")
+ except IOError:
+ model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
+ return
+ except IOError:
+ model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
+ return
+ ### find idea for derivate in obprocess.computeGradV
+ ### TODO: read file + create varfs like txtmodel
+ ###Â we have \Delta H(u),
+ ### H_comps, Dp_{ai}H_comps
+ self.varfs=[None]*(self.nbComps*(self.dS.study.nbParam+1))
+ self.varfsTm1=[None]*self.nbComps
+ #for the comp number numComp \in {0,..,self.nbComps-1}:
+ #self.varfs[numComp*(self.nbComps+1)] contains the varf F_numComp for the rhs
+ #self.varfs[numComp*(self.nbComps+1)+(j+1)] contains the derivative d_jF_numComp, j \in {0,..,self.nbComps-1}
+ strkeys=r"[\s\+\-\*/\=\(\)\^]+"
+ for numComps in range(self.nbComps):
+ ##READ SECOND ORDER TERM
+ strline=f.readline().split('=')[1].strip()
+ strline=strForDolf(strline, self.unknownName)
+ strline=strline.replace('^','**')
+ ##replace param
+ for count_p in range(self.dS.paramD):
+ strline=strline.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
+ model.myPrint(" read txt model F"+txtModel.unknownName[numComps]+"="+strline)
+ #### replace COVARIABLES
+ strline=updadeVarfWithCovariables(strline,self.dS.covariables)
+ model.myPrint(strline)
+ #### BUILD VARF USING FK
+ code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+"*self.dS.v*self.dS.dx)"
+ exec(compile(code, 'sumstring', 'exec'))
+ #### BUILD VARF USING UTM1
+ code="self.varfsTm1["+str(numComps)+"]="+self._subFXForm(strline,"utm1")+"*self.dS.v*self.dS.dx)"
+ model.myPrint(code)
+ exec(compile(code, 'sumstring', 'exec'))
+ ## READ DERIVATIVES PARAM
+ for countP in range(self.dS.paramD):
+ strline=fd.readline().split('=')[1].strip()
+ strline=strForDolf(strline,self.unknownName)
+ #strline=strline.replace('^','**')
+ ##replace param
+ for tmpCountP in enumerate(self.dS.paramD):
+ strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountp)+']')
+ strline=updadeVarfWithCovariables(strline,self.dS.covariables)
+ model.myPrint("txtModel: read txt model d"+txtModel.unknownName[k]+"F"+txtModel.unknownName[numComps]+"="+strline)
+ code="self.varfs["+str(numComps)+"*(self.dS.study.nbParam+1)+"+str(k)+"+1]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
+ exec(compile(code, 'sumstring', 'exec'))
def updateBlockMat(self,compi,compj):
+<<<<<<< HEAD
if (compi==compj):
return self.block
def _computeBlockAdjoint(self):
@@ -173,6 +261,9 @@ class txtDiffusionModel(model):
self.bilinearForm=form(aux)
self.blockAdjoint=matrix(assemble_matrix(self.bilinearForm))
>>>>>>> 19312c4 (quick update, not finish)
+=======
+ return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)+compj+1]))
+>>>>>>> ef8e663 (update init of txtDiffusionModel, not finish)
def updateBlockMatAdjoint(self,compi,compj):
if (compi==compj):
return self.blockAdjoint
--
GitLab
From 62e184df8c406b0570b86fe341f9cedf34490948 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Fri, 31 May 2024 16:27:54 +0200
Subject: [PATCH 16/23] remove all message from gitlab
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 125 --------------------
1 file changed, 125 deletions(-)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index b816abc..4daaa63 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -9,26 +9,8 @@ from mse.linAlgTools import matrix,vector
class diffusionModel(model):
def __init__(self,aDs,aDiffCoef=1) -> None:
super().__init__(aDs)
-<<<<<<< HEAD
-<<<<<<< HEAD
-<<<<<<< HEAD
aDs.addSecondOrderTerm(self)
self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
-=======
- self.coefName=aDiffCoef
->>>>>>> 57c81bb (Deletes modif in diffuisionModel that added a parameter. NOTE: this modif will be moved to another class)
-=======
-=======
->>>>>>> ef8e663 (update init of txtDiffusionModel, not finish)
- self.coefName=aDiffCoef
-=======
- ### NOTE Test: we define diffusionModel with param
-=======
->>>>>>> 5c06eb8 (update init of txtDiffusionModel, not finish)
- aDs.addSecondOrderTerm(self)
- self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
->>>>>>> 3f748d9 (quick update, not finish)
->>>>>>> 19312c4 (quick update, not finish)
self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
self.block=None
self.blockAdjoint=None
@@ -72,7 +54,6 @@ class txtDiffusionModel(model):
self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
self._computeBlock()
self.updateLinearisedMat()
-<<<<<<< HEAD
def _computeBlock(self):
aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
self.bilinearForm=form(aux)
@@ -115,7 +96,6 @@ class txtDiffusionModel(model):
except IOError:
model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
return
-<<<<<<< HEAD
### find idea for derivate in obprocess.computeGradV
### TODO: read file + create varfs like txtmodel
###Â we have \Delta H(u),
@@ -159,111 +139,6 @@ class txtDiffusionModel(model):
exec(compile(code, 'sumstring', 'exec'))
def updateBlockMat(self,compi,compj):
return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)+compj+1]))
-=======
- except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
- return
- ### find idea for derivate in obprocess.computeGradV
- ### TODO: read file + create varfs like txtmodel
- def truc():
- pass
-=======
->>>>>>> ef8e663 (update init of txtDiffusionModel, not finish)
- def _computeBlock(self):
- aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
- self.bilinearForm=form(aux)
- self.block=matrix(assemble_matrix(self.bilinearForm))
- # we check if file exist
- if (exists(aTxtFile+".txt")):
- self.txtFile=aTxtFile
- else:
- self.txtFile=computerEnv.modelsDir+aTxtFile #aTxtFile form is "affineGrowth"
- model.myPrint("Looking for model "+self.txtFile)
- try:
- f = open(self.txtFile+".txt", "r")
- except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+".txt")
- return
- # we check if the derivate of parameter file exist
- try:
- fd = open(self.txtFile+"DVP.txt", "r")
- except FileNotFoundError:
- # if didn't exist, we create it
- self.dvFileExist=False
- print("Création fichier "+self.txtFile+"DVP.tx")
- fd = open(self.txtFile+"DVP.txt", "w") #open file in mod write
- for numComps in range(self.nbComps):
- strF=f.readline()
- strF=strF.split('=')[1].strip()#.replace("^","**")
- for p in range(self.aDs.study.nbParam):
- grad1s=diff(strF,'p'+str(p+1))
- grad1=str(grad1s).replace("**","^")
- #write in fic
- fd.write("dp"+str(p+1)+"f"+txtModel.unknownName[numComps]+"="+grad1+"\n")
- f.close()
- fd.close()
- f = open(self.txtFile+".txt", "r")
- try:
- fd = open(self.txtFile+"DVP.txt", "r")
- except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
- return
- except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
- return
- ### find idea for derivate in obprocess.computeGradV
- ### TODO: read file + create varfs like txtmodel
- ###Â we have \Delta H(u),
- ### H_comps, Dp_{ai}H_comps
- self.varfs=[None]*(self.nbComps*(self.dS.study.nbParam+1))
- self.varfsTm1=[None]*self.nbComps
- #for the comp number numComp \in {0,..,self.nbComps-1}:
- #self.varfs[numComp*(self.nbComps+1)] contains the varf F_numComp for the rhs
- #self.varfs[numComp*(self.nbComps+1)+(j+1)] contains the derivative d_jF_numComp, j \in {0,..,self.nbComps-1}
- strkeys=r"[\s\+\-\*/\=\(\)\^]+"
- for numComps in range(self.nbComps):
- ##READ SECOND ORDER TERM
- strline=f.readline().split('=')[1].strip()
- strline=strForDolf(strline, self.unknownName)
- strline=strline.replace('^','**')
- ##replace param
- for count_p in range(self.dS.paramD):
- strline=strline.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
- model.myPrint(" read txt model F"+txtModel.unknownName[numComps]+"="+strline)
- #### replace COVARIABLES
- strline=updadeVarfWithCovariables(strline,self.dS.covariables)
- model.myPrint(strline)
- #### BUILD VARF USING FK
- code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+"*self.dS.v*self.dS.dx)"
- exec(compile(code, 'sumstring', 'exec'))
- #### BUILD VARF USING UTM1
- code="self.varfsTm1["+str(numComps)+"]="+self._subFXForm(strline,"utm1")+"*self.dS.v*self.dS.dx)"
- model.myPrint(code)
- exec(compile(code, 'sumstring', 'exec'))
- ## READ DERIVATIVES PARAM
- for countP in range(self.dS.paramD):
- strline=fd.readline().split('=')[1].strip()
- strline=strForDolf(strline,self.unknownName)
- #strline=strline.replace('^','**')
- ##replace param
- for tmpCountP in enumerate(self.dS.paramD):
- strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountp)+']')
- strline=updadeVarfWithCovariables(strline,self.dS.covariables)
- model.myPrint("txtModel: read txt model d"+txtModel.unknownName[k]+"F"+txtModel.unknownName[numComps]+"="+strline)
- code="self.varfs["+str(numComps)+"*(self.dS.study.nbParam+1)+"+str(k)+"+1]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
- exec(compile(code, 'sumstring', 'exec'))
- def updateBlockMat(self,compi,compj):
-<<<<<<< HEAD
- if (compi==compj):
- return self.block
- def _computeBlockAdjoint(self):
- aux=dot(grad(self.dS.u),grad(self.diffCoef*self.dS.v))*self.dS.dx
- self.bilinearForm=form(aux)
- self.blockAdjoint=matrix(assemble_matrix(self.bilinearForm))
->>>>>>> 19312c4 (quick update, not finish)
-=======
- return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)+compj+1]))
->>>>>>> ef8e663 (update init of txtDiffusionModel, not finish)
def updateBlockMatAdjoint(self,compi,compj):
if (compi==compj):
return self.blockAdjoint
--
GitLab
From e5631201534f5d56d0ad7a254b743766b5300fef Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Wed, 5 Jun 2024 13:51:03 +0200
Subject: [PATCH 17/23] Patch txtDispersionModel: the 2ndOrderTerm is update
with COVARIABLE
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 92 ++++++++++++---------
trunk/src/mse/obsProcess.py | 2 +-
2 files changed, 52 insertions(+), 42 deletions(-)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 4daaa63..6adc05e 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -1,4 +1,4 @@
-from mse.MODELS.models import model, strForDolf
+from mse.MODELS.models import model, strForDolf, updadeVarfWithCovariables
import mse.dynamicalSystem
from dolfinx.fem import Constant,form
from petsc4py.PETSc import ScalarType
@@ -10,7 +10,6 @@ class diffusionModel(model):
def __init__(self,aDs,aDiffCoef=1) -> None:
super().__init__(aDs)
aDs.addSecondOrderTerm(self)
- self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
self.block=None
self.blockAdjoint=None
@@ -40,11 +39,10 @@ class txtDiffusionModel(model):
""" Class used to create a second order term with a txt file
"""
+ unknownName="abcdefgh"
def __init__(self,aDs, aTxtFile) -> None:
super().__init__(aDs)
aDs.addSecondOrderTerm(self)
- self.coefName=aDiffCoef#FIXME: rm all call to this then rm this line
- self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
self.block=None
self.blockAdjoint=None
self.bilinearForm=None
@@ -77,14 +75,18 @@ class txtDiffusionModel(model):
self.dvFileExist=False
print("Création fichier "+self.txtFile+"DVP.tx")
fd = open(self.txtFile+"DVP.txt", "w") #open file in mod write
+ # loop for each comp
for numComps in range(self.nbComps):
- strF=f.readline()
- strF=strF.split('=')[1].strip()#.replace("^","**")
- for p in range(self.aDs.study.nbParam):
- grad1s=diff(strF,'p'+str(p+1))
- grad1=str(grad1s).replace("**","^")
- #write in fic
- fd.write("dp"+str(p+1)+"f"+txtModel.unknownName[numComps]+"="+grad1+"\n")
+ #loop for each dimension(=2)
+ for strDim in ['x','y']:
+ strF=f.readline()
+ strF=strF.split('=')[1].strip()#.replace("^","**")
+ #loop for each param
+ for p in range(self.aDs.study.nbParam):
+ grad1s=diff(strF,'p'+str(p+1))
+ grad1=str(grad1s).replace("**","^")
+ #write in fic
+ fd.write("dp"+strDim+str(p+1)+"f"+txtModel.unknownName[numComps]+"="+grad1+"\n")
f.close()
fd.close()
f = open(self.txtFile+".txt", "r")
@@ -98,47 +100,55 @@ class txtDiffusionModel(model):
return
### find idea for derivate in obprocess.computeGradV
### TODO: read file + create varfs like txtmodel
- ###Â we have \Delta H(u),
- ### H_comps, Dp_{ai}H_comps
- self.varfs=[None]*(self.nbComps*(self.dS.study.nbParam+1))
- self.varfsTm1=[None]*self.nbComps
+ nbDim = 2
+ self.varfs=[None]*(self.nbComps*nbDim*(self.dS.study.nbParam+1))
+ #self.varfsTm1=[None]*self.nbComps
#for the comp number numComp \in {0,..,self.nbComps-1}:
#self.varfs[numComp*(self.nbComps+1)] contains the varf F_numComp for the rhs
#self.varfs[numComp*(self.nbComps+1)+(j+1)] contains the derivative d_jF_numComp, j \in {0,..,self.nbComps-1}
strkeys=r"[\s\+\-\*/\=\(\)\^]+"
for numComps in range(self.nbComps):
##READ SECOND ORDER TERM
- strline=f.readline().split('=')[1].strip()
- strline=strForDolf(strline, self.unknownName)
- strline=strline.replace('^','**')
- ##replace param
- for count_p in range(self.dS.paramD):
- strline=strline.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
- model.myPrint(" read txt model F"+txtModel.unknownName[numComps]+"="+strline)
- #### replace COVARIABLES
- strline=updadeVarfWithCovariables(strline,self.dS.covariables)
- model.myPrint(strline)
- #### BUILD VARF USING FK
- code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+"*self.dS.v*self.dS.dx)"
- exec(compile(code, 'sumstring', 'exec'))
- #### BUILD VARF USING UTM1
- code="self.varfsTm1["+str(numComps)+"]="+self._subFXForm(strline,"utm1")+"*self.dS.v*self.dS.dx)"
- model.myPrint(code)
- exec(compile(code, 'sumstring', 'exec'))
- ## READ DERIVATIVES PARAM
- for countP in range(self.dS.paramD):
- strline=fd.readline().split('=')[1].strip()
- strline=strForDolf(strline,self.unknownName)
- #strline=strline.replace('^','**')
+ # loop for each dimension (2 dim)
+ for numDim in range(nbDim)
+ strline=f.readline().split('=')[1].strip()
+ strline=strForDolf(strline, self.unknownName)
+ strline=strline.replace('^','**')
##replace param
- for tmpCountP in enumerate(self.dS.paramD):
- strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountp)+']')
+ for count_p in range(self.dS.paramD):
+ strline=strline.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
+ #### replace COVARIABLES
strline=updadeVarfWithCovariables(strline,self.dS.covariables)
- model.myPrint("txtModel: read txt model d"+txtModel.unknownName[k]+"F"+txtModel.unknownName[numComps]+"="+strline)
- code="self.varfs["+str(numComps)+"*(self.dS.study.nbParam+1)+"+str(k)+"+1]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
+ model.myPrint(strline)
+ #### BUILD VARF USING FK
+ code="self.varfs["+str(numComps)+"*"+str(numDim+1)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+"*self.dS.v*self.dS.dx)"
+ print("num index in varf =",numComps*(numDim+1)*(self.nbComps+1))
exec(compile(code, 'sumstring', 'exec'))
+ ##### BUILD VARF USING UTM1
+ #code="self.varfsTm1["+str(numComps)+"]="+self._subFXForm(strline,"utm1")+"*self.dS.v*self.dS.dx)"
+ #model.myPrint(code)
+ #exec(compile(code, 'sumstring', 'exec'))
+ ## READ DERIVATIVES PARAM
+ # loop for each param
+ for countP in range(self.dS.paramD):
+ strline=fd.readline().split('=')[1].strip()
+ strline=strForDolf(strline,self.unknownName)
+ ##replace param
+ for tmpCountP in enumerate(self.dS.paramD):
+ strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountp)+']')
+ #replace COVARIABLES
+ strline=updadeVarfWithCovariables(strline,self.dS.covariables)
+ model.myPrint("txtModel: read txt model d"+txtModel.unknownName[k]+"F"+txtModel.unknownName[numComps]+"="+strline)
+ code="self.varfs["+str(numComps)+"*"+str(numDim+1)+"*(self.dS.study.nbParam+1)+"+str(k)+"+1]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
+ print("num index in derivate param =",numComps*(numDim+1)*(self.dS.study.nbParam+1)+k+1)
+ exec(compile(code, 'sumstring', 'exec'))
def updateBlockMat(self,compi,compj):
return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)+compj+1]))
def updateBlockMatAdjoint(self,compi,compj):
if (compi==compj):
return self.blockAdjoint
+ #TODO: surcharge function "getlineraized" to add update
+ def getLinearisedMat(self):
+ super().updateLinearisedMat()
+ super().getLinearisedMat()
+ # def get(): super().update, super.get()
diff --git a/trunk/src/mse/obsProcess.py b/trunk/src/mse/obsProcess.py
index 5b3e8dd..0a0d341 100644
--- a/trunk/src/mse/obsProcess.py
+++ b/trunk/src/mse/obsProcess.py
@@ -199,7 +199,7 @@ class obsProcess:
#read u_a
code="self.daF.x.array[:]="+'+'.join(strdFp[i])
exec(compile(code,'sumstring', 'exec'))
- curF[i] = assemble_scalar(form(comp*self.daF*self.dS.dx))
+ curF[i] = assemble_scalar(form(comp*self.daF*self.dS.dx)) #+ derive du terme de 2nd ordre par rapport a p TODO
### End V1
### V2
# NOTE: cette version demande de creer un vecteur. L'impementation est à l'abandon
--
GitLab
From 9b7c38d3982947b97418ab7d1531b8e3e22aaded Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Wed, 5 Jun 2024 13:52:16 +0200
Subject: [PATCH 18/23] add test for new class txtDiffModel
---
trunk/tests/test_secndOrdrTermTxt.py | 87 ++++++++++++++++++++++++++++
1 file changed, 87 insertions(+)
create mode 100644 trunk/tests/test_secndOrdrTermTxt.py
diff --git a/trunk/tests/test_secndOrdrTermTxt.py b/trunk/tests/test_secndOrdrTermTxt.py
new file mode 100644
index 0000000..3e2f7d3
--- /dev/null
+++ b/trunk/tests/test_secndOrdrTermTxt.py
@@ -0,0 +1,87 @@
+
+from mse.study import study
+from mse.dynamicalSystem import dynamicalSystem,computeMass
+from mse.DISPERSIONS.diffusionModel import diffusionModel
+from mse.MODELS.models import txtModel
+from mse.system import system
+from mse.simulator import simulator
+from mse.timeDoOneStep import timeDoOneStep
+from mse.TIMESCHEMES.timeScheme import implicitLinearTimeScheme,implicitNonLinearTimeScheme
+from mse.toolsEnv import computerEnv,MSEPrint
+from mse.explorTrajectory import explorTrajectory
+
+import numpy as np
+import pytest
+from os import mkdir
+from os.path import isdir
+
+
+
+### NOTE: without sourceTerm ? only 2nd order term read in file ?
+##Â test with constant sourceTerm + other.
+
+
+#un modele
+nameModele="preyPred"
+nameSendOrdTerm="sOTConst"
+sourceModele = computerEnv.modelsDir+nameModele+".txt"
+nbComps=2
+nbparam=2
+eps=1e-6
+eps2=1e-2
+
+ #une etude
+studName="TEST_SNDORDTERM"
+aStudy=study(studName,nbparam)
+aStudy.setGeometry("SQUARE",0.4)
+ #un systeme
+aSystem=system(aStudy)
+ #un systeme dynamique
+ds2d1=dynamicalSystem(aStudy,nbComps,"surface1.xdmf",nbComps=nbComps)
+
+m2=txtModel(ds2d1,nameModele)
+ #m2.computeMatOnce=1
+ #m2.computeRhsOnce=0
+
+aTS=None
+#schema en temps implicite lineaire
+aTS=implicitLinearTimeScheme()
+aTS.computeMatOnce=1
+
+simulator(0,1,aStudy)
+aTStep=timeDoOneStep(aStudy,aTS,0.1)
+mass=np.zeros((ds2d1.nbComps),dtype=np.double)
+
+#add diffusion motion
+aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
+m=txtModel(ds2d1,nameModel)
+#
+#aTS=None
+## select the time algorithm
+#if (testNLTS):
+# aTS=implicitNonLinearTimeScheme()
+#else:
+# aTS=implicitLinearTimeScheme()
+##simu in t \in [0,1]
+#simulator(0,1,aStudy)
+##based on fix time step 0.05
+#aTStep=timeDoOneStep(aStudy,aTS,0.05)
+##run the simulation
+#aStudy.simulator.doSimu()
+#
+#
+##for visu, export simulation steps to vtk
+#aStudy.simulator.exportSimuToVtk()
+#
+##loop on simu to compute mass
+#import numpy as np
+#mass=np.zeros((ds2d1.nbComps),dtype=np.double)
+#aExplorTraj=explorTrajectory(aStudy,aStudy.simuDir)
+#aExplorTraj.rewind()
+#while(aExplorTraj.replay()):
+# computeMass(ds2d1,mass)
+# print("mass "+str(aExplorTraj.curStep)+" = ",end="")
+# for m in mass:
+# print(m,end=", ")
+# print("")
+#aStudy.end()
--
GitLab
From 2b73b91903a3f8c78a0effc992503a06b442d056 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Fri, 7 Jun 2024 15:39:22 +0200
Subject: [PATCH 19/23] rm void line in test
---
trunk/tests/test_StudyParam2.py | 4 ----
1 file changed, 4 deletions(-)
diff --git a/trunk/tests/test_StudyParam2.py b/trunk/tests/test_StudyParam2.py
index 5bedb02..944658a 100644
--- a/trunk/tests/test_StudyParam2.py
+++ b/trunk/tests/test_StudyParam2.py
@@ -80,7 +80,3 @@ def test_Param2(result, index, param, trueParam):
print("final mass = ",mass[0])
assert r
-
-
-
-
--
GitLab
From 8d658e16d4558439b1199aa5123da699d68b2383 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Fri, 7 Jun 2024 15:44:52 +0200
Subject: [PATCH 20/23] Patch class txtDiffModel: work with parameter
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 241 +++++++++++++-------
1 file changed, 153 insertions(+), 88 deletions(-)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 6adc05e..78eee11 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -1,10 +1,14 @@
-from mse.MODELS.models import model, strForDolf, updadeVarfWithCovariables
+from mse.MODELS.models import model, strForDolf, formatCovariables
import mse.dynamicalSystem
-from dolfinx.fem import Constant,form
+from dolfinx.fem import Constant,form,assemble_scalar
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import assemble_matrix
-from ufl import grad,dot
+from ufl import grad,dot, inner
from mse.linAlgTools import matrix,vector
+from mse.toolsEnv import computerEnv, MSEPrint, MSEPrint_0
+from os.path import exists
+from sympy import diff
+import numpy as np
class diffusionModel(model):
def __init__(self,aDs,aDiffCoef=1) -> None:
@@ -23,6 +27,8 @@ class diffusionModel(model):
def _computeBlock(self):
aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
self.bilinearForm=form(aux)
+ print("aux = ",aux)
+ print("bilineaire form in diff model(cst): ",self.bilinearForm)
self.block=matrix(assemble_matrix(self.bilinearForm))
def updateBlockMat(self,compi,compj):
if (compi==compj):
@@ -40,9 +46,17 @@ class txtDiffusionModel(model):
"""
unknownName="abcdefgh"
+ def _subFXForm(self,strModel,X):
+ for j in range(self.nbComps):
+ strModel=strModel.replace("$"+txtDiffusionModel.unknownName[j],"self.dS."+X+".comps["+str(j)+"]")
+ code="form(("+strModel+")"
+ model.myPrint("X ->"+code)
+ return code
+
def __init__(self,aDs, aTxtFile) -> None:
super().__init__(aDs)
aDs.addSecondOrderTerm(self)
+ self.txtFile=aTxtFile
self.block=None
self.blockAdjoint=None
self.bilinearForm=None
@@ -53,102 +67,153 @@ class txtDiffusionModel(model):
self._computeBlock()
self.updateLinearisedMat()
def _computeBlock(self):
- aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
- self.bilinearForm=form(aux)
- self.block=matrix(assemble_matrix(self.bilinearForm))
# we check if file exist
- if (exists(aTxtFile+".txt")):
- self.txtFile=aTxtFile
- else:
- self.txtFile=computerEnv.modelsDir+aTxtFile #aTxtFile form is "affineGrowth"
- model.myPrint("Looking for model "+self.txtFile)
- try:
- f = open(self.txtFile+".txt", "r")
- except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+".txt")
- return
- # we check if the derivate of parameter file exist
+ if (exists(self.txtFile+".txt")):
+ self.txtFile=self.txtFile
+ else:
+ self.txtFile=computerEnv.modelsDir+self.txtFile #self.txtFile form is "affineGrowth"
+ model.myPrint("Looking for model "+self.txtFile)
+ try:
+ f = open(self.txtFile+".txt", "r")
+ except IOError:
+ model.myPrint("IOError in txtDiffModel when reading "+self.txtFile+".txt")
+ return
+ # we check if the derivate of parameter file exist
+ try:
+ fd = open(self.txtFile+"DVP.txt", "r")
+ except FileNotFoundError:
+ # if didn't exist, we create it
+ self.dvFileExist=False
+ print("Création fichier "+self.txtFile+"DVP.tx")
+ fd = open(self.txtFile+"DVP.txt", "w") #open file in mod write
+ # loop for each comp
+ for numComps in range(self.nbComps):
+ #loop for each dimension(=2)
+ for strDim in ['x','y']:
+ strF=f.readline()
+ strF=strF.split('=')[1].strip()#.replace("^","**")
+ #loop for each param
+ for p in range(self.dS.study.nbParam):
+ grad1s=diff(strF,'p'+str(p+1))
+ grad1=str(grad1s).replace("**","^")
+ #write in fic
+ fd.write("dp"+strDim+str(p+1)+"f"+txtDiffusionModel.unknownName[numComps]+"="+grad1+"\n")
+ f.close()
+ fd.close()
+ f = open(self.txtFile+".txt", "r")
try:
fd = open(self.txtFile+"DVP.txt", "r")
- except FileNotFoundError:
- # if didn't exist, we create it
- self.dvFileExist=False
- print("Création fichier "+self.txtFile+"DVP.tx")
- fd = open(self.txtFile+"DVP.txt", "w") #open file in mod write
- # loop for each comp
- for numComps in range(self.nbComps):
- #loop for each dimension(=2)
- for strDim in ['x','y']:
- strF=f.readline()
- strF=strF.split('=')[1].strip()#.replace("^","**")
- #loop for each param
- for p in range(self.aDs.study.nbParam):
- grad1s=diff(strF,'p'+str(p+1))
- grad1=str(grad1s).replace("**","^")
- #write in fic
- fd.write("dp"+strDim+str(p+1)+"f"+txtModel.unknownName[numComps]+"="+grad1+"\n")
- f.close()
- fd.close()
- f = open(self.txtFile+".txt", "r")
- try:
- fd = open(self.txtFile+"DVP.txt", "r")
- except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
- return
except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+"DVP.txt")
+ model.myPrint("IOError in txtDiffusionModel when reading "+self.txtFile+"DVP.txt")
return
- ### find idea for derivate in obprocess.computeGradV
- ### TODO: read file + create varfs like txtmodel
- nbDim = 2
- self.varfs=[None]*(self.nbComps*nbDim*(self.dS.study.nbParam+1))
- #self.varfsTm1=[None]*self.nbComps
- #for the comp number numComp \in {0,..,self.nbComps-1}:
- #self.varfs[numComp*(self.nbComps+1)] contains the varf F_numComp for the rhs
- #self.varfs[numComp*(self.nbComps+1)+(j+1)] contains the derivative d_jF_numComp, j \in {0,..,self.nbComps-1}
- strkeys=r"[\s\+\-\*/\=\(\)\^]+"
- for numComps in range(self.nbComps):
- ##READ SECOND ORDER TERM
- # loop for each dimension (2 dim)
- for numDim in range(nbDim)
- strline=f.readline().split('=')[1].strip()
- strline=strForDolf(strline, self.unknownName)
- strline=strline.replace('^','**')
+ except IOError:
+ model.myPrint("IOError in txtDiffusionModel when reading "+self.txtFile+"DVP.txt")
+ return
+ ### find idea for derivate in obprocess.computeGradV
+ ### TODO: read file + create varfs like txtmodel
+ nbDim = 2
+ self.varfs=[None]*(self.nbComps*(nbDim*self.dS.study.nbParam+1))
+ #self.varfsTm1=[None]*self.nbComps
+ #for the comp number numComp \in {0,..,self.nbComps-1}:
+ #self.varfs[numComp*(self.nbComps+1)] contains the varf F_numComp for the rhs
+ #self.varfs[numComp*(self.nbComps+1)+(j+1)] contains the derivative d_jF_numComp, j \in {0,..,self.nbComps-1}
+ strkeys=r"[\s\+\-\*/\=\(\)\^]+"
+ for numComps in range(self.nbComps):
+ ##READ SECOND ORDER TERM
+ # Dxx
+ strlineX=f.readline().split('=')[1].strip()
+ strlineX=strForDolf(strlineX, self.unknownName)
+ strlineX=strlineX.replace('^','**')
+ # Dyy
+ strlineY=f.readline().split('=')[1].strip()
+ strlineY=strForDolf(strlineY, self.unknownName)
+ strlineY=strlineY.replace('^','**')
+ ##replace param
+ for count_p in range(len(self.dS.paramD)):
+ strlineX=strlineX.replace('p'+str(count_p+1),'self.dS.study.param['+str(count_p)+']')
+ strlineY=strlineY.replace('p'+str(count_p+1),'self.dS.study.param['+str(count_p)+']')
+ #### replace COVARIABLES
+ strlineX=formatCovariables(strlineX,"self.dS.dicCov[\"", "\"]")
+ strlineY=formatCovariables(strlineY,"self.dS.dicCov[\"", "\"]")
+ #### BUILD VARF USING FK
+ # v1: with Grad[i]
+ strline = strlineX + "*grad(self.dS.u)[0]*grad(self.dS.v)[0]*self.dS.dx+"+strlineY+"*grad(self.dS.u)[1]*grad(self.dS.v)[1]*self.dS.dx"
+ model.myPrint(strline)
+ print("varfs = ",strline)
+ code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+")"
+ print("num index in varf =",numComps*(self.nbComps+1))
+ print("code = ",code)
+ exec(compile(code, 'sumstring', 'exec'))
+ ##### BUILD VARF USING UTM1
+ #code="self.varfsTm1["+str(numComps)+"]="+self._subFXForm(strline,"utm1")+"*self.dS.v*self.dS.dx)"
+ #model.myPrint(code)
+ #exec(compile(code, 'sumstring', 'exec'))
+ # v2: with D= const(np.array[[Dxx,0],[0,Dyy]])
+ #Dxy=[None]*2
+ #code="Dxy[0]="+strlineX
+ #print("code Dxx = ",code)
+ #exec(compile(code, 'sumstring', 'exec'))
+ #code="Dxy[1]="+strlineY
+ #print("code Dyy = ",code)
+ #exec(compile(code, 'sumstring', 'exec'))
+ #print("Dxx = "+str(Dxy[0])+", Dyy = ",Dxy[1])
+ ### Mat
+ ##D = Constant(self.dS.dolfinMesh, np.array([[Dxy[0], 0], [0, Dxy[1]]], dtype=np.float64))
+ ###Â Vec
+ #D = Constant(self.dS.dolfinMesh, np.array([Dxy[0], Dxy[1]], dtype=np.float64))
+ ### R
+ ##D = Constant(self.dS.dolfinMesh, self.dS.study.param[0])
+ #strline = "dot(D*grad(self.dS.u)+ self.dS.u*grad(D),grad(self.dS.v))*self.dS.dx"
+ #code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+")"
+ #print("num index in varf =",numComps*(self.nbComps+1))
+ #print("code = ",code)
+ ##exec(compile(code, 'sumstring', 'exec'))
+ #self.varfs[0*(self.nbComps+1)]=form((dot(dot(D,grad(self.dS.u))+ dot(self.dS.u,grad(D)),grad(self.dS.v))*self.dS.dx))
+ ## READ DERIVATIVES PARAM
+ # loop for each dim (=2)
+ for numDim in range(nbDim):
+ # loop for each param
+ for countP in range(len(self.dS.paramD)):
+ strline=fd.readline().split('=')[1].strip()
+ strline=strForDolf(strline,self.unknownName)
##replace param
- for count_p in range(self.dS.paramD):
- strline=strline.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
- #### replace COVARIABLES
- strline=updadeVarfWithCovariables(strline,self.dS.covariables)
- model.myPrint(strline)
- #### BUILD VARF USING FK
- code="self.varfs["+str(numComps)+"*"+str(numDim+1)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+"*self.dS.v*self.dS.dx)"
- print("num index in varf =",numComps*(numDim+1)*(self.nbComps+1))
+ for tmpCountP,tmpP in enumerate(self.dS.paramD):
+ strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountP)+']')
+ #replace COVARIABLES
+ strline=formatCovariables(strline,"self.dS.dicCov[\"", "\"]")
+ model.myPrint("txtDiffusionModel: read txt model d"+txtDiffusionModel.unknownName[numComps]+"F"+txtDiffusionModel.unknownName[numComps]+"="+strline)
+ code="self.varfs["+str(numComps)+"*(self.dS.study.nbParam+1)+"+str(numComps)+"+1+"+str(numDim)+"]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
+ print("codeDvp = ",code)
+ print("num index in derivate param =",numComps*(self.dS.study.nbParam+1)+numComps+1+(numDim))
exec(compile(code, 'sumstring', 'exec'))
- ##### BUILD VARF USING UTM1
- #code="self.varfsTm1["+str(numComps)+"]="+self._subFXForm(strline,"utm1")+"*self.dS.v*self.dS.dx)"
- #model.myPrint(code)
- #exec(compile(code, 'sumstring', 'exec'))
- ## READ DERIVATIVES PARAM
- # loop for each param
- for countP in range(self.dS.paramD):
- strline=fd.readline().split('=')[1].strip()
- strline=strForDolf(strline,self.unknownName)
- ##replace param
- for tmpCountP in enumerate(self.dS.paramD):
- strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountp)+']')
- #replace COVARIABLES
- strline=updadeVarfWithCovariables(strline,self.dS.covariables)
- model.myPrint("txtModel: read txt model d"+txtModel.unknownName[k]+"F"+txtModel.unknownName[numComps]+"="+strline)
- code="self.varfs["+str(numComps)+"*"+str(numDim+1)+"*(self.dS.study.nbParam+1)+"+str(k)+"+1]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
- print("num index in derivate param =",numComps*(numDim+1)*(self.dS.study.nbParam+1)+k+1)
- exec(compile(code, 'sumstring', 'exec'))
+ #### create block,NOTE: useless
+ #Dx = self.varfs[0]
+ #Dy = self.varfs[1]
+ #print("we will print each comp of aux:")
+ #print("Dx: ",type(Dx))
+ #print("Dx= ",Dx)
+ #print("grad(self.dS.u)[0]: ",type(grad(self.dS.u)[0]))
+ #print("grad(self.dS.u)[0]= ",grad(self.dS.u)[0])
+ #print(",grad(self.dS.v)[0] ",type(grad(self.dS.v)[0]))
+ #print("self.dS.dx, ",type(self.dS.dx))
+ #print("Dy, ",type(Dy))
+ #print("grad(self.dS.u)[1], ",type(grad(self.dS.u)[1]))
+ #print("grad(self.dS.v)[1], ",type(grad(self.dS.v)[1]))
+ #print("self.dS.dx, ",type(self.dS.dx))
+ #matD = Constant(self.dS.dolfinMesh, np.array([ [Dx,0], [0,Dy] ], dtype=np.float64))
+ #self.bilinearForm = Dx*grad(self.dS.u)[0]*grad(self.dS.v)[0]*self.dS.dx+Dy*grad(self.dS.u)[1]*grad(self.dS.v)[1]*self.dS.dx
+ #self.bilinearForm=form(aux)
+ #self.block=matrix(assemble_matrix(self.bilinearForm))
def updateBlockMat(self,compi,compj):
- return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)+compj+1]))
+ if (compi==compj):
+ #return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)]))
+ return matrix(assemble_matrix(self.varfs[0]))
def updateBlockMatAdjoint(self,compi,compj):
if (compi==compj):
return self.blockAdjoint
- #TODO: surcharge function "getlineraized" to add update
+ ### surcharge function "getlineraized" to add update
def getLinearisedMat(self):
super().updateLinearisedMat()
- super().getLinearisedMat()
+ return self.mat
+ #NOTE: super()getLinearizedMat didn't work....
# def get(): super().update, super.get()
--
GitLab
From 3d50ab7db4d2046b64d2f8abea9f954825567a2b Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Fri, 7 Jun 2024 15:45:28 +0200
Subject: [PATCH 21/23] Test for class txtDiffModel
---
trunk/tests/test_secndOrdrTermTxt.py | 83 ++++++++++++++--------------
1 file changed, 43 insertions(+), 40 deletions(-)
diff --git a/trunk/tests/test_secndOrdrTermTxt.py b/trunk/tests/test_secndOrdrTermTxt.py
index 3e2f7d3..b4cb2b8 100644
--- a/trunk/tests/test_secndOrdrTermTxt.py
+++ b/trunk/tests/test_secndOrdrTermTxt.py
@@ -1,7 +1,6 @@
-
from mse.study import study
from mse.dynamicalSystem import dynamicalSystem,computeMass
-from mse.DISPERSIONS.diffusionModel import diffusionModel
+from mse.DISPERSIONS.diffusionModel import diffusionModel, txtDiffusionModel
from mse.MODELS.models import txtModel
from mse.system import system
from mse.simulator import simulator
@@ -22,13 +21,15 @@ from os.path import isdir
#un modele
-nameModele="preyPred"
+nameModele="affineGrowth"
nameSendOrdTerm="sOTConst"
sourceModele = computerEnv.modelsDir+nameModele+".txt"
-nbComps=2
+nbComps=1
+WITH_AFFINE_GROWTH=1
nbparam=2
eps=1e-6
eps2=1e-2
+testNLTS=0
#une etude
studName="TEST_SNDORDTERM"
@@ -37,12 +38,9 @@ aStudy.setGeometry("SQUARE",0.4)
#un systeme
aSystem=system(aStudy)
#un systeme dynamique
-ds2d1=dynamicalSystem(aStudy,nbComps,"surface1.xdmf",nbComps=nbComps)
+ds2d1=dynamicalSystem(aStudy,2,"surface1.xdmf",nbComps=nbComps)
-m2=txtModel(ds2d1,nameModele)
- #m2.computeMatOnce=1
- #m2.computeRhsOnce=0
-
+aStudy.setParam([0.,1.])
aTS=None
#schema en temps implicite lineaire
aTS=implicitLinearTimeScheme()
@@ -54,34 +52,39 @@ mass=np.zeros((ds2d1.nbComps),dtype=np.double)
#add diffusion motion
aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
-m=txtModel(ds2d1,nameModel)
-#
-#aTS=None
-## select the time algorithm
-#if (testNLTS):
-# aTS=implicitNonLinearTimeScheme()
-#else:
-# aTS=implicitLinearTimeScheme()
-##simu in t \in [0,1]
-#simulator(0,1,aStudy)
-##based on fix time step 0.05
-#aTStep=timeDoOneStep(aStudy,aTS,0.05)
-##run the simulation
-#aStudy.simulator.doSimu()
-#
-#
-##for visu, export simulation steps to vtk
-#aStudy.simulator.exportSimuToVtk()
-#
-##loop on simu to compute mass
-#import numpy as np
-#mass=np.zeros((ds2d1.nbComps),dtype=np.double)
-#aExplorTraj=explorTrajectory(aStudy,aStudy.simuDir)
-#aExplorTraj.rewind()
-#while(aExplorTraj.replay()):
-# computeMass(ds2d1,mass)
-# print("mass "+str(aExplorTraj.curStep)+" = ",end="")
-# for m in mass:
-# print(m,end=", ")
-# print("")
-#aStudy.end()
+#aDiffModel=diffusionModel(ds2d1)
+#aDiffModel.setDiffCoef(1.)
+if (WITH_AFFINE_GROWTH):
+ m=txtModel(ds2d1,nameModele)
+ m.computeMatOnce=1
+ m.computeRhsOnce=1
+
+aTS=None
+# select the time algorithm
+if (testNLTS):
+ aTS=implicitNonLinearTimeScheme()
+else:
+ aTS=implicitLinearTimeScheme()
+#simu in t \in [0,1]
+simulator(0,1,aStudy)
+#based on fix time step 0.05
+aTStep=timeDoOneStep(aStudy,aTS,0.05)
+#run the simulation
+aStudy.simulator.doSimu()
+
+print("param = ",aStudy.param)
+#for visu, export simulation steps to vtk
+aStudy.simulator.exportSimuToVtk()
+
+#loop on simu to compute mass
+import numpy as np
+mass=np.zeros((ds2d1.nbComps),dtype=np.double)
+aExplorTraj=explorTrajectory(aStudy,aStudy.simuDir)
+aExplorTraj.rewind()
+while(aExplorTraj.replay()):
+ computeMass(ds2d1,mass)
+ print("mass "+str(aExplorTraj.curStep)+" = ",end="")
+ for m in mass:
+ print(m,end=", ")
+ print("")
+aStudy.end()
--
GitLab
From a6d56e1924e02ff236006a40965a4aa97f2a6627 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Mon, 10 Jun 2024 14:33:25 +0200
Subject: [PATCH 22/23] Add features to compute Jacobian if 2ndOrdTerm contain
parameter, to test
---
trunk/src/mse/DISPERSIONS/diffusionModel.py | 49 +++++++--------------
trunk/src/mse/obsProcess.py | 18 ++++++--
2 files changed, 31 insertions(+), 36 deletions(-)
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 78eee11..7a70258 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -5,7 +5,7 @@ from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import assemble_matrix
from ufl import grad,dot, inner
from mse.linAlgTools import matrix,vector
-from mse.toolsEnv import computerEnv, MSEPrint, MSEPrint_0
+from mse.toolsEnv import computerEnv, MSEPrint, MSEPrint_0
from os.path import exists
from sympy import diff
import numpy as np
@@ -141,7 +141,7 @@ class txtDiffusionModel(model):
model.myPrint(strline)
print("varfs = ",strline)
code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+")"
- print("num index in varf =",numComps*(self.nbComps+1))
+ print("num index in varf =",numComps*(self.nbComps+1))
print("code = ",code)
exec(compile(code, 'sumstring', 'exec'))
##### BUILD VARF USING UTM1
@@ -165,7 +165,7 @@ class txtDiffusionModel(model):
##D = Constant(self.dS.dolfinMesh, self.dS.study.param[0])
#strline = "dot(D*grad(self.dS.u)+ self.dS.u*grad(D),grad(self.dS.v))*self.dS.dx"
#code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+")"
- #print("num index in varf =",numComps*(self.nbComps+1))
+ #print("num index in varf =",numComps*(self.nbComps+1))
#print("code = ",code)
##exec(compile(code, 'sumstring', 'exec'))
#self.varfs[0*(self.nbComps+1)]=form((dot(dot(D,grad(self.dS.u))+ dot(self.dS.u,grad(D)),grad(self.dS.v))*self.dS.dx))
@@ -175,35 +175,20 @@ class txtDiffusionModel(model):
# loop for each param
for countP in range(len(self.dS.paramD)):
strline=fd.readline().split('=')[1].strip()
- strline=strForDolf(strline,self.unknownName)
- ##replace param
- for tmpCountP,tmpP in enumerate(self.dS.paramD):
- strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountP)+']')
- #replace COVARIABLES
- strline=formatCovariables(strline,"self.dS.dicCov[\"", "\"]")
- model.myPrint("txtDiffusionModel: read txt model d"+txtDiffusionModel.unknownName[numComps]+"F"+txtDiffusionModel.unknownName[numComps]+"="+strline)
- code="self.varfs["+str(numComps)+"*(self.dS.study.nbParam+1)+"+str(numComps)+"+1+"+str(numDim)+"]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
- print("codeDvp = ",code)
- print("num index in derivate param =",numComps*(self.dS.study.nbParam+1)+numComps+1+(numDim))
- exec(compile(code, 'sumstring', 'exec'))
- #### create block,NOTE: useless
- #Dx = self.varfs[0]
- #Dy = self.varfs[1]
- #print("we will print each comp of aux:")
- #print("Dx: ",type(Dx))
- #print("Dx= ",Dx)
- #print("grad(self.dS.u)[0]: ",type(grad(self.dS.u)[0]))
- #print("grad(self.dS.u)[0]= ",grad(self.dS.u)[0])
- #print(",grad(self.dS.v)[0] ",type(grad(self.dS.v)[0]))
- #print("self.dS.dx, ",type(self.dS.dx))
- #print("Dy, ",type(Dy))
- #print("grad(self.dS.u)[1], ",type(grad(self.dS.u)[1]))
- #print("grad(self.dS.v)[1], ",type(grad(self.dS.v)[1]))
- #print("self.dS.dx, ",type(self.dS.dx))
- #matD = Constant(self.dS.dolfinMesh, np.array([ [Dx,0], [0,Dy] ], dtype=np.float64))
- #self.bilinearForm = Dx*grad(self.dS.u)[0]*grad(self.dS.v)[0]*self.dS.dx+Dy*grad(self.dS.u)[1]*grad(self.dS.v)[1]*self.dS.dx
- #self.bilinearForm=form(aux)
- #self.block=matrix(assemble_matrix(self.bilinearForm))
+ if strline == '0':
+ pass
+ else:
+ strline=strForDolf(strline,self.unknownName)
+ ##replace param
+ for tmpCountP,tmpP in enumerate(self.dS.paramD):
+ strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountP)+']')
+ #replace COVARIABLES
+ strline=formatCovariables(strline,"self.dS.dicCov[\"", "\"]")
+ model.myPrint("txtDiffusionModel: read txt model d"+txtDiffusionModel.unknownName[numComps]+"F"+txtDiffusionModel.unknownName[numComps]+"="+strline)
+ code="self.varfs["+str(numComps)+"*(self.dS.study.nbParam+1)+"+str(numComps)+"+1+"+str(numDim)+"]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
+ print("codeDvp = ",code)
+ print("num index in derivate param =",numComps*(self.dS.study.nbParam+1)+numComps+1+(numDim))
+ exec(compile(code, 'sumstring', 'exec'))
def updateBlockMat(self,compi,compj):
if (compi==compj):
#return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)]))
diff --git a/trunk/src/mse/obsProcess.py b/trunk/src/mse/obsProcess.py
index 0a0d341..7c449b1 100644
--- a/trunk/src/mse/obsProcess.py
+++ b/trunk/src/mse/obsProcess.py
@@ -14,6 +14,7 @@ from mse.systemAdjoint import systemAdjoint
from mse.TIMESCHEMES.timeScheme import implicitLinearBackwardTimeScheme
from mse.timeDoOneStep import timeDoOneStepBackward
from mse.simulator import simulatorBackward
+from mse.toolsEnv import MSEPrint, MSEPrint_0
from os.path import exists
from os import mkdir
class obsProcess:
@@ -23,6 +24,8 @@ class obsProcess:
"""
def __init__(self, aDS, aDirURef=None)->None:
#Dynamical system
+ self.self.myPrint=MSEPrint_0
+ #self.self.myPrint=MSEPrint
self.dS=aDS
self.dS.setObsProcess(self)
self.index = self.dS.indexInSystem
@@ -78,7 +81,7 @@ class obsProcess:
try:
f = open(sourceTerm.txtFile+".txt", "r")
except IOError:
- model.myPrint("IOError in txtModel when reading "+self.txtFile+".txt")
+ self.myPrint("IOError in txtModel when reading "+self.txtFile+".txt")
return
try:
fp = open(sourceTerm.txtFile+"Dp"+".txt", "r")
@@ -96,7 +99,7 @@ class obsProcess:
#write in fic
fp.write("d"+"p"+str(k+1)+"f"+sourceTerm.unknownName[numComp]+"="+grad1+"\n")
except IOError:
- model.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
+ self.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
return
f.close()
fp.close()
@@ -117,7 +120,7 @@ class obsProcess:
# try:
# fp = open(sourceTerm.txtFile+"Dp"+".txt", "r")
# except IOError:
-# model.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
+# self.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
# return
# for numComp in range(self.dS.nbComps):
# for numParam in range(self.dS.study.nbParam):
@@ -148,7 +151,7 @@ class obsProcess:
try:
fp = open(sourceTerm.txtFile+"Dp"+".txt", "r")
except IOError:
- model.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
+ self.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
return
for numComp in range(self.dS.nbComps):
for numParam in range(self.dS.study.nbParam):
@@ -173,6 +176,7 @@ class obsProcess:
pastT=self.dS.study.simulatorBackward.t0
pastF=np.zeros(self.dS.study.nbParam)
curF =np.zeros(self.dS.study.nbParam)
+ ### todo: READ SECONDoRDERtER.varfs[??], if != None add in grad
# TODO: reprendre comme model: tableau de ufl, 1 dim de taille nbparam*???
while(l_aExplorTrajAdjoint.replay()):
#adjoint is load in ds.utm1 after call of l_aExplorTrajAdjoint.replay
@@ -194,11 +198,17 @@ class obsProcess:
# loop for each param
# comp contain P, u_a is in trajState.comps[]
#tmp_count+=1#FIXME
+ tmpDPSOT = 0
for i in range(self.dS.study.nbParam):
###Â V1
#read u_a
code="self.daF.x.array[:]="+'+'.join(strdFp[i])
exec(compile(code,'sumstring', 'exec'))
+ # derivate of 2ndOrderTerm
+ if(self.dS.secondOrderTerm.varfs[1+i+count] = None):tmpDPSOT = 0
+ else:
+ self.myPrint("count = ",1+i+count)
+ tmpDPSOT = assemble_scalar(form(-self.dS.secondOrderTerm.varfs[1+i+count]*inner(grad(self.daF),grad(comp)*self.dS.dx)))
curF[i] = assemble_scalar(form(comp*self.daF*self.dS.dx)) #+ derive du terme de 2nd ordre par rapport a p TODO
### End V1
### V2
--
GitLab
From 2faa99d33310279ba5c366791709be0f70a134b3 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Wed, 12 Jun 2024 11:04:36 +0200
Subject: [PATCH 23/23] patch txtDiffModel+obsProcess to compute gradV with
parameter + test
---
trunk/BD_MODELS/sOTConst.txt | 2 ++
trunk/BD_MODELS/sOTConst2.txt | 2 ++
trunk/src/mse/DISPERSIONS/diffusionModel.py | 22 +++++++-----
trunk/src/mse/obsProcess.py | 38 +++++++++++++++++----
trunk/tests/test_backwardModel.py | 2 +-
5 files changed, 50 insertions(+), 16 deletions(-)
create mode 100644 trunk/BD_MODELS/sOTConst.txt
create mode 100644 trunk/BD_MODELS/sOTConst2.txt
diff --git a/trunk/BD_MODELS/sOTConst.txt b/trunk/BD_MODELS/sOTConst.txt
new file mode 100644
index 0000000..a98da38
--- /dev/null
+++ b/trunk/BD_MODELS/sOTConst.txt
@@ -0,0 +1,2 @@
+fax=1*p1
+fay=p2
diff --git a/trunk/BD_MODELS/sOTConst2.txt b/trunk/BD_MODELS/sOTConst2.txt
new file mode 100644
index 0000000..42071d1
--- /dev/null
+++ b/trunk/BD_MODELS/sOTConst2.txt
@@ -0,0 +1,2 @@
+fax=p3
+fay=p3
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 7a70258..ef13327 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -112,7 +112,7 @@ class txtDiffusionModel(model):
### find idea for derivate in obprocess.computeGradV
### TODO: read file + create varfs like txtmodel
nbDim = 2
- self.varfs=[None]*(self.nbComps*(nbDim*self.dS.study.nbParam+1))
+ self.varfs=[None]*(self.nbComps*(nbDim*self.dS.study.nbParam+1))#FIXME
#self.varfsTm1=[None]*self.nbComps
#for the comp number numComp \in {0,..,self.nbComps-1}:
#self.varfs[numComp*(self.nbComps+1)] contains the varf F_numComp for the rhs
@@ -130,8 +130,8 @@ class txtDiffusionModel(model):
strlineY=strlineY.replace('^','**')
##replace param
for count_p in range(len(self.dS.paramD)):
- strlineX=strlineX.replace('p'+str(count_p+1),'self.dS.study.param['+str(count_p)+']')
- strlineY=strlineY.replace('p'+str(count_p+1),'self.dS.study.param['+str(count_p)+']')
+ strlineX=strlineX.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
+ strlineY=strlineY.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
#### replace COVARIABLES
strlineX=formatCovariables(strlineX,"self.dS.dicCov[\"", "\"]")
strlineY=formatCovariables(strlineY,"self.dS.dicCov[\"", "\"]")
@@ -140,8 +140,8 @@ class txtDiffusionModel(model):
strline = strlineX + "*grad(self.dS.u)[0]*grad(self.dS.v)[0]*self.dS.dx+"+strlineY+"*grad(self.dS.u)[1]*grad(self.dS.v)[1]*self.dS.dx"
model.myPrint(strline)
print("varfs = ",strline)
- code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+")"
- print("num index in varf =",numComps*(self.nbComps+1))
+ code="self.varfs["+str(numComps)+"*(2*self.dS.study.nbParam+1)]="+self._subFXForm(strline,"uk")+")"
+ print("num index in varf =",numComps*(2*self.dS.study.nbParam+1))
print("code = ",code)
exec(compile(code, 'sumstring', 'exec'))
##### BUILD VARF USING UTM1
@@ -169,26 +169,32 @@ class txtDiffusionModel(model):
#print("code = ",code)
##exec(compile(code, 'sumstring', 'exec'))
#self.varfs[0*(self.nbComps+1)]=form((dot(dot(D,grad(self.dS.u))+ dot(self.dS.u,grad(D)),grad(self.dS.v))*self.dS.dx))
+ #
## READ DERIVATIVES PARAM
+ # NOTE: the value didn't depend of 'u', so didn't need 'self.subxForm'
+ # NOTE: we will use 'lambda' to create function ref to param.
# loop for each dim (=2)
for numDim in range(nbDim):
# loop for each param
for countP in range(len(self.dS.paramD)):
strline=fd.readline().split('=')[1].strip()
+ print("num index for derivate param =",numComps*(2*self.dS.study.nbParam+1)+countP+1+(numDim)*(self.dS.study.nbParam))
+ print("strline==\"0\":",strline=="0")
if strline == '0':
pass
else:
strline=strForDolf(strline,self.unknownName)
##replace param
for tmpCountP,tmpP in enumerate(self.dS.paramD):
- strline=strline.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountP)+']')
+ strline=strline.replace('p'+str(tmpCountP+1),'self.dS.study.param['+str(tmpCountP)+']')
#replace COVARIABLES
+ #tmpDPSOT = assemble_scalar(form((tmpDPSOTX + tmpDPSOTY)))
strline=formatCovariables(strline,"self.dS.dicCov[\"", "\"]")
model.myPrint("txtDiffusionModel: read txt model d"+txtDiffusionModel.unknownName[numComps]+"F"+txtDiffusionModel.unknownName[numComps]+"="+strline)
- code="self.varfs["+str(numComps)+"*(self.dS.study.nbParam+1)+"+str(numComps)+"+1+"+str(numDim)+"]="+self._subFXForm(strline,"uk")+"*self.dS.u*self.dS.v*self.dS.dx)"
+ code="self.varfs["+str(numComps)+"*(2*self.dS.study.nbParam+1)+"+str(tmpCountP)+"+1+"+str(numDim)+"*(self.dS.study.nbParam)]=lambda: "+strline
print("codeDvp = ",code)
- print("num index in derivate param =",numComps*(self.dS.study.nbParam+1)+numComps+1+(numDim))
exec(compile(code, 'sumstring', 'exec'))
+ for i in self.varfs:print (i)
def updateBlockMat(self,compi,compj):
if (compi==compj):
#return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)]))
diff --git a/trunk/src/mse/obsProcess.py b/trunk/src/mse/obsProcess.py
index 7c449b1..622e605 100644
--- a/trunk/src/mse/obsProcess.py
+++ b/trunk/src/mse/obsProcess.py
@@ -24,8 +24,8 @@ class obsProcess:
"""
def __init__(self, aDS, aDirURef=None)->None:
#Dynamical system
- self.self.myPrint=MSEPrint_0
- #self.self.myPrint=MSEPrint
+ #self.myPrint=MSEPrint_0
+ self.myPrint=MSEPrint
self.dS=aDS
self.dS.setObsProcess(self)
self.index = self.dS.indexInSystem
@@ -178,6 +178,8 @@ class obsProcess:
curF =np.zeros(self.dS.study.nbParam)
### todo: READ SECONDoRDERtER.varfs[??], if != None add in grad
# TODO: reprendre comme model: tableau de ufl, 1 dim de taille nbparam*???
+ tmpDPSOTX = 0
+ tmpDPSOTY = 0
while(l_aExplorTrajAdjoint.replay()):
#adjoint is load in ds.utm1 after call of l_aExplorTrajAdjoint.replay
# coef si integration rectangle ou trapeze
@@ -198,18 +200,40 @@ class obsProcess:
# loop for each param
# comp contain P, u_a is in trajState.comps[]
#tmp_count+=1#FIXME
- tmpDPSOT = 0
for i in range(self.dS.study.nbParam):
###Â V1
#read u_a
code="self.daF.x.array[:]="+'+'.join(strdFp[i])
exec(compile(code,'sumstring', 'exec'))
+ # 2 DIM
# derivate of 2ndOrderTerm
- if(self.dS.secondOrderTerm.varfs[1+i+count] = None):tmpDPSOT = 0
+ # X
+ if (not hasattr(self.dS.secondOrderTerm,"varfs")):tmpDPSOT=0
else:
- self.myPrint("count = ",1+i+count)
- tmpDPSOT = assemble_scalar(form(-self.dS.secondOrderTerm.varfs[1+i+count]*inner(grad(self.daF),grad(comp)*self.dS.dx)))
- curF[i] = assemble_scalar(form(comp*self.daF*self.dS.dx)) #+ derive du terme de 2nd ordre par rapport a p TODO
+ indX = count*(2*self.dS.study.nbParam+1)+i+1
+ indY = count*(2*self.dS.study.nbParam+1)+i+1+self.dS.study.nbParam
+ DpDx = self.dS.secondOrderTerm.varfs[indX]
+ DpDy = self.dS.secondOrderTerm.varfs[indY]
+ if(DpDx == None):tmpDPSOTX = 0*self.dS.dx
+ else:
+ tmpDPSOTX = DpDx()*grad(self.dS.trajState.comps[count])[0]*grad(comp)[0]*self.dS.dx
+ # Y
+ if(DpDy == None):tmpDPSOTY = 0*self.dS.dx
+ else:
+ tmpDPSOTY = DpDy()*grad(self.dS.trajState.comps[count])[1]*grad(comp)[1]*self.dS.dx
+ #print("type of tmpDPSOTX: ",type(tmpDPSOTX))
+ #print("(tmpDPSOTX): ",tmpDPSOTX)
+ #print("type of form(tmpDPSOTX): ",type(form(tmpDPSOTX)))
+ #print("type of :tmpDPSOTY ",type(tmpDPSOTY))
+ #print("(tmpDPSOTY): ",tmpDPSOTY)
+ #print("type of form(tmpDPSOTY): ",type(form(tmpDPSOTY)))
+ #print("type of :form(comp*self.daF*self.dS.dx) ",type(form(comp*self.daF*self.dS.dx)))
+ #print("type of : form((tmpDPSOTX + tmpDPSOTY))) ",type(form((tmpDPSOTX + tmpDPSOTY))))
+ #print("type of : assemble_scalar(form((tmpDPSOTX + tmpDPSOTY)))) ",type(assemble_scalar(form((tmpDPSOTX + tmpDPSOTY)))))
+ #print("assemble_scalar(form((tmpDPSOTX + tmpDPSOTY)))) ",assemble_scalar(form((tmpDPSOTX + tmpDPSOTY))))
+ tmpDPSOT = assemble_scalar(form(tmpDPSOTX + tmpDPSOTY))
+ print("tmpDPSOT: ",tmpDPSOT)
+ curF[i] = assemble_scalar(form(comp*self.daF*self.dS.dx)) + tmpDPSOT
### End V1
### V2
# NOTE: cette version demande de creer un vecteur. L'impementation est à l'abandon
diff --git a/trunk/tests/test_backwardModel.py b/trunk/tests/test_backwardModel.py
index 9e6af36..511c2df 100644
--- a/trunk/tests/test_backwardModel.py
+++ b/trunk/tests/test_backwardModel.py
@@ -90,7 +90,7 @@ quaDiff=quadraticDiff(ds2d1)
# We test with different epsilon and parameters.
print("Param a_ref =",aStudy.getParam)
-@pytest.mark.parametrize("defaultParam",[([0.5,0.1]),([1.5,1.])])
+@pytest.mark.parametrize("defaultParam",[([0.5,0.1]),([1.5,1.]),([1.2,.6])])
def test_GradWithFiniteDifference(defaultParam):
"""
test the approximation of GradV mse with finite Difference.
--
GitLab