diff --git a/trunk/BD_MODELS/sOTConst.txt b/trunk/BD_MODELS/sOTConst.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a98da38d3263464c3c3b0a838c1588159138c15f
--- /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 0000000000000000000000000000000000000000..42071d10f6545365b3637689a3416bbf6b760f27
--- /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 8184738722cf35efeca7b36b38e8bfa46ee61d77..ef13327c5885de651bd9797d4e898b2b6b017d7d 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -1,15 +1,19 @@
-from mse.MODELS.models import model
+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:
super().__init__(aDs)
- self.coefName=aDiffCoef
+ aDs.addSecondOrderTerm(self)
self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
self.block=None
self.blockAdjoint=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):
@@ -34,3 +40,171 @@ class diffusionModel(model):
def updateBlockMatAdjoint(self,compi,compj):
if (compi==compj):
return self.blockAdjoint
+
+class txtDiffusionModel(model):
+ """ Class used to create a second order term with a txt file
+
+ """
+ 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
+ self._computeBlock()
+ self.updateLinearisedMat()
+ def setDiffCoef(self,aTxtFile):
+ self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
+ self._computeBlock()
+ self.updateLinearisedMat()
+ def _computeBlock(self):
+ # we check if 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 IOError:
+ model.myPrint("IOError in txtDiffusionModel when reading "+self.txtFile+"DVP.txt")
+ return
+ 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))#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
+ #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.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[\"", "\"]")
+ #### 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)+"*(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
+ #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
+ # 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.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)+"*(2*self.dS.study.nbParam+1)+"+str(tmpCountP)+"+1+"+str(numDim)+"*(self.dS.study.nbParam)]=lambda: "+strline
+ print("codeDvp = ",code)
+ 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)]))
+ return matrix(assemble_matrix(self.varfs[0]))
+ def updateBlockMatAdjoint(self,compi,compj):
+ if (compi==compj):
+ return self.blockAdjoint
+ ### surcharge function "getlineraized" to add update
+ def getLinearisedMat(self):
+ super().updateLinearisedMat()
+ return self.mat
+ #NOTE: super()getLinearizedMat didn't work....
+ # def get(): super().update, super.get()
diff --git a/trunk/src/mse/obsProcess.py b/trunk/src/mse/obsProcess.py
index 5b3e8dd1cfb43077e5709130354c484e9c972343..622e60506efd116ee39d403244541763ef3a4c61 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.myPrint=MSEPrint_0
+ 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,7 +176,10 @@ 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*???
+ 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
@@ -199,7 +205,35 @@ 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))
+ # 2 DIM
+ # derivate of 2ndOrderTerm
+ # X
+ if (not hasattr(self.dS.secondOrderTerm,"varfs")):tmpDPSOT=0
+ else:
+ 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_StudyParam2.py b/trunk/tests/test_StudyParam2.py
index 5bedb024d3da5e5f92628116e045683d680f65a1..944658aa47a8fe3cd04e7975395589cf475ecd40 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
-
-
-
-
diff --git a/trunk/tests/test_backwardModel.py b/trunk/tests/test_backwardModel.py
index 9e6af3619f847a54faeafd70c8d4795e78c5c7f0..511c2df3a41bd3a976228f78a306eb34be813fb7 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.
diff --git a/trunk/tests/test_secndOrdrTermTxt.py b/trunk/tests/test_secndOrdrTermTxt.py
new file mode 100644
index 0000000000000000000000000000000000000000..b4cb2b81783a80c78cb861d5ff470d5bdd8dd813
--- /dev/null
+++ b/trunk/tests/test_secndOrdrTermTxt.py
@@ -0,0 +1,90 @@
+from mse.study import study
+from mse.dynamicalSystem import dynamicalSystem,computeMass
+from mse.DISPERSIONS.diffusionModel import diffusionModel, txtDiffusionModel
+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="affineGrowth"
+nameSendOrdTerm="sOTConst"
+sourceModele = computerEnv.modelsDir+nameModele+".txt"
+nbComps=1
+WITH_AFFINE_GROWTH=1
+nbparam=2
+eps=1e-6
+eps2=1e-2
+testNLTS=0
+
+ #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,2,"surface1.xdmf",nbComps=nbComps)
+
+aStudy.setParam([0.,1.])
+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)
+#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()