diff --git a/trunk/BD_MODELS/sOTConst3.txt b/trunk/BD_MODELS/sOTConst3.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c47554185466dbaa9c229e0a120500d774f200a8
--- /dev/null
+++ b/trunk/BD_MODELS/sOTConst3.txt
@@ -0,0 +1,2 @@
+fax=1
+fay=p1
diff --git a/trunk/BD_MODELS/sOTConst4.txt b/trunk/BD_MODELS/sOTConst4.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e3ba67b97edd85c476cea3d53dc6e262fdccf17b
--- /dev/null
+++ b/trunk/BD_MODELS/sOTConst4.txt
@@ -0,0 +1,2 @@
+fax=p2
+fay=1
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 7ac089d6f438e1a2271dd0deee446e6172784309..b9bdab6be1d5b516eef99877e13508098e7beb72 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -29,7 +29,7 @@ class diffusionModel(model):
self.updateLinearisedMat()
self.nbDim = 2
self.varfs=[None for i in range(self.nbComps*(self.dS.study.nbParam+1))]
-
+
def setDiffCoef(self,aDiffCoef):
self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
self._computeBlock()
@@ -55,7 +55,7 @@ class diffusionModel(model):
class txtDiffusionModel(model):
""" Class used to create a second order term with a txt file
-
+
NOTE: only avaible for compment in 2d
:param dynamicalSystem aDs: class allows applying a PDE to a geometry
:param str strModel: file contain diffModel.The txt file must be this form:
@@ -72,6 +72,19 @@ class txtDiffusionModel(model):
model.myPrint("X ->"+code)
return code
+ def _execStr(self, strEquation, index):
+ """Function will take a string wich contain equation, and do execution to transform it in form
+ """
+ ##replace param
+ for tmpCountP,tmpP in enumerate(self.dS.paramD):
+ strEquation=strEquation.replace('p'+str(tmpCountP+1),'self.dS.paramD['+str(tmpCountP)+']')
+ #replace COVARIABLES
+ strEquation=formatCovariables(strEquation,"self.dS.dicCov[\"", "\"]")
+ #### BUILD VARF
+ code="self.varfs["+str(index)+"]=form("+strEquation+")"
+ exec(compile(code, 'sumstring', 'exec'))
+
+
def __init__(self,aDs, aTxtFile) -> None:
super().__init__(aDs)
aDs.addSecondOrderTerm(self)
@@ -167,10 +180,8 @@ class txtDiffusionModel(model):
# 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)
+ index = numComps*(2*self.dS.study.nbParam+1)
+ code="self.varfs["+ str(index) +"]="+self._subFXForm(strline,"uk")+")"
exec(compile(code, 'sumstring', 'exec'))
## READ DERIVATIVES PARAM
# NOTE: we will save form, if 0, we save NULL or 0
@@ -190,37 +201,26 @@ class txtDiffusionModel(model):
strlineY=strForDolf(strlineY,self.unknownName)
# index
index = numComps*(2*self.dS.study.nbParam+1)+countP+1
- print("index for derivate: ",index)
- print("strline x = "+strlineX+"\nstrliney = "+strlineY)
##Â Test if deivate is Null
if strlineX == "0":
if strlineY =="0":
## if both der ==Null: we replace Varfs[index] by None
self.varfs[index]=None
- else:
+ else:
#Â Only dim x ==0
- strline ="-grad("+strlineY+"*self.dS.trajState.comps["+str(numComps)+"])[1]*grad(self.dS.bufAdjoint.comps["+str(numComps)+"])[1])*self.dS.dx"
+ strline ="-grad("+strlineY+"*self.dS.trajState.comps["+str(numComps)+"])[1]*grad(self.dS.bufAdjoint.comps["+str(numComps)+"])[1]*self.dS.dx"
+ # Function to transform str in form and put it in varfs
+ self._execStr(strline, index)
elif strlineY == "0":
#Â Only dim y ==0
strline = "-grad("+strlineX+"*self.dS.trajState.comps["+str(numComps)+"])[0]*grad(self.dS.bufAdjoint.comps["+str(numComps)+"])[0]*self.dS.dx"
+ # Function to transform str in form and put it in varfs
+ self._execStr(strline, index)
else:
# both != 0
strline = "(-grad("+strlineX+"*self.dS.trajState.comps["+str(numComps)+"])[0]*grad(self.dS.bufAdjoint.comps["+str(numComps)+"])[0] - grad("+strlineY+"*self.dS.trajState.comps["+str(numComps)+"])[1]*grad(self.dS.bufAdjoint.comps["+str(numComps)+"])[1])*self.dS.dx"
- ##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[\"", "\"]")
- #### BUILD VARF USING
- #tmpDPSOT = assemble_scalar(form((tmpDPSOTX + tmpDPSOTY)))
- #model.myPrint("txtDiffusionModel: read txt model d"+txtDiffusionModel.unknownName[numComps]+"F"+txtDiffusionModel.unknownName[numComps]+"="+strlineX + " + "+strlineY)
- ### NOTE use _subXForm
- # Assemble#TODO: _subXform for form
- code="self.varfs["+str(index)+"]=form("+strline+")"
- #code="self.varfs["+str(numComps)+"*(2*self.dS.study.nbParam+1)+"+str(tmpCountP)+"+1*(self.dS.study.nbParam)]="+strline
- print("codeDvp = ",code)
- exec(compile(code, 'sumstring', 'exec'))
- for i in self.varfs:print (i)
+ # Function to transform str in form and put it in varfs
+ self._execStr(strline, index)
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 ef018317c26b114e59547d34fa900533b4e5514d..868366394359344bb9fe2dc7a78bf472e9e0b436 100644
--- a/trunk/src/mse/obsProcess.py
+++ b/trunk/src/mse/obsProcess.py
@@ -24,6 +24,10 @@ class obsProcess:
This class contain the definition of main function for an obsprocess.
"""
def _subFXForm(self,strModel,X):
+ """function use for exec of str in Code.
+
+ we transform a string of Model in a string wich contain form and the dX
+ """
#NOTE: move this function ? Used in txtModel, txtDiffModel
unknownName="abcdefgh"
for j in range(self.dS.nbComps):
@@ -62,8 +66,9 @@ class obsProcess:
self.simuBackward = simulatorBackward(t1,t0,self.dS.study)
self.tDoStepBackward = timeDoOneStepBackward(self.dS.study, self.timeSchemeBackward)
self.timeSchemeBackward.theta = 1
- # Tab wich contain derivate
- self.varDerParam =[0*self.dS.dx for i in range (self.dS.study.nbParam)] #[None]*(self.dS.study.nbParam)
+ # Tab wich contain the sum of all derivate of ST with respect to param (size of nbParam)
+ #this is initiate in computeGradV()
+ self.varDerParam = None
def computeV(self):
""" Method to compute the observation process for the current study parameters.
"""
@@ -73,7 +78,7 @@ class obsProcess:
"""
pass
def computeGradV(self):
- """By default, it builds the adjoint state via systemAdjoint class and build the grad of the obsprocess using <P.\partial_aF(a,ua)e_i>.
+ """By default, it builds the adjoint state via systemAdjoint class and build the grad of the obsprocess using <P.\\partial_aF(a,ua)e_i>.
read file from adjoint simu directory and create gradV.
"""
@@ -89,69 +94,48 @@ class obsProcess:
for sourceTerm in self.dS.sourcesTerm:
sourceTerm.enableAdjoint()
self.dS.secondOrderTerm.enableAdjoint()
+ ##3
+ ## Check if varDerParam is init, else create it
+ if self.varDerParam != None:pass
+ else:
+ self.varDerParam = [None]*(self.dS.study.nbParam)
+ # tab of str we wil excecute
+ strdFp = [ '' for _ in range(self.dS.study.nbParam)] #tab while contain str of dFp for each comps, =daF
+ # loop for each ST
+ for sourceTerm in self.dS.sourcesTerm:
+ #tmp str, value of '' or '+'
+ # We read derivate of ST
+ try:
+ fp = open(sourceTerm.txtFile+"Dp"+".txt", "r")
+ except IOError:
+ 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):
+ # read the line
+ l_linedFp = fp.readline().split('=')[1].strip()
+ #Test if the derivate is !=0
+ if l_linedFp != "0":
+ # replace compment for dolfinix
+ l_linedFp = strForDolf(l_linedFp,sourceTerm.unknownName)
+ # replace param
+ for count,p in enumerate(self.dS.study.param):
+ l_linedFp=l_linedFp.replace('p'+str(count+1),'self.dS.paramD['+str(count)+']')
+ # Test if strFD is empty or not => if had to add '+' or not
+ if (strdFp[numParam] == ''):strdFp[numParam]=l_linedFp
+ else:strdFp[numParam]='+'.join(strdFp[numParam],l_linedFp)
+ else:pass
+ fp.close()
+ #EXCECUTION
+ for count,strline in enumerate(strdFp):
+ if strline =='':pass
+ else:
+ code = "self.varDerParam ["+str(count)+"]="+self._subFXForm(strline,"trajState")+"*self.dS.bufAdjoint.comps["+str(numComp)+"]*self.dS.dx)"
+ print("in obsProcess: code = "+code)
+ exec(compile(code,'sumstring', 'exec'))
##2: create vardp := \partial_aF(a,ua)
- # V1: On garde en memoires les strlines puis remplacer "u" par "ua". Problemes si beaucoup de parametres.
- # V2: using same method than model with varfs: create empty list size of nbSourceTerm*nbParam*nbcomps, then
- # fill in the table the val of integral.
-# ##### BEGIN V2
- #NOTE: cette version demande le passage d'un var ufl.form.Form en dolfinx.fem.function.Function.
- # Abandon pour le moment
-# # create the vec vardfa:
-# nbST = len(self.dS.sourcesTerm)
-# self.vardfa = [None]*(nbST*self.dS.nbComps*self.dS.study.nbParam)
-# print("size of vardfa = ",nbST*self.dS.nbComps*self.dS.study.nbParam)
-# # loop for each sourceTerm of ds
-# numST = 0
-# for sourceTerm in self.dS.sourcesTerm:
-# try:
-# fp = open(sourceTerm.txtFile+"Dp"+".txt", "r")
-# except IOError:
-# 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):
-# strline_dFp = fp.readline().split('=')[1].strip()
-# strline_dFp = strForDolf(strline_dFp,sourceTerm.unknownName)
-# for count,p in enumerate(self.dS.study.param):
-# strline_dFp =strline_dFp.replace('p'+str(count+1),'self.dS.study.param['+str(count)+']')
-# for countComp in range(self.dS.nbComps):
-# strline_dFp = strline_dFp.replace("$"+sourceTerm.unknownName[countComp],"self.dS.trajState.comps["+str(countComp)+"].x.array[:]")
-# tmpIndST = numST*self.dS.nbComps*self.dS.study.nbParam
-# tmpIndComp = numComp*self.dS.study.nbParam
-# tmpIndparam = numParam
-# tmpInd = tmpIndST+tmpIndComp+tmpIndparam
-# print("tmpIndST =",tmpIndST)
-# print("tmpIndComp =",tmpIndComp)
-# print("tmpIndparam =",tmpIndparam)
-# code = "self.vardfa["+str(tmpInd)+"]="+strline_dFp+"*self.dS.u*self.dS.v*self.dS.dx"
-# print("str of vardfa = ",code)
-# exec(compile(code, 'sumstring', 'exec'))
-# #print("vardfa = ",self.vardfa[tmpInd])
-# numST+=1
-# fp.close()
-# ###### END
###V1: with str
- strdFp = [] #tab while contain str of dFp for each comps, =daF
- for _ in range(self.dS.study.nbParam):strdFp.append([])
- for sourceTerm in self.dS.sourcesTerm:
- try:
- fp = open(sourceTerm.txtFile+"Dp"+".txt", "r")
- except IOError:
- 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):
- l_linedFp = fp.readline().split('=')[1].strip()
- l_linedFp = strForDolf(l_linedFp,sourceTerm.unknownName)
- for count,p in enumerate(self.dS.study.param):
- l_linedFp=l_linedFp.replace('p'+str(count+1),'self.dS.study.param['+str(count)+']')
- #NOTE: delete 2 pro lines pour utiliser subXForm ?
- for countComp in range(self.dS.nbComps):
- l_linedFp = l_linedFp.replace("$"+sourceTerm.unknownName[countComp],"self.dS.trajState.comps["+str(countComp)+"].x.array[:]")
- #print("l_linedFp = "+l_linedFp)
- strdFp[numParam].append(l_linedFp)
- fp.close()
- #NOTE: l_strdFp is 2d list, l_strdFp[i][j] return the partial derivate of F_j model of a_i param
+ #NOTE: l_strdFp is 2d list, l_strdFp[i][j] return the partial derivate with respect of F_j model of a_i param
# to use it: dFa_i = "+".join(strdFp[i])pastT
l_aExplorTrajAdjoint=explorTrajectory(self.dS.study,self.dS.study.adjointDir)
l_aExplorTrajAdjoint.rewind()
@@ -163,7 +147,6 @@ class obsProcess:
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()):
@@ -177,27 +160,29 @@ class obsProcess:
#Update ua
self.explor_u.interpolToTime(curT)
l_dtAdj = abs(curT - pastT)
- # loop for each sourceterm
- for numST in range(len(self.dS.sourcesTerm)):
- # loop for each comps
- for count,comp in enumerate(self.dS.bufAdjoint.comps):
- #NOTE: on suppose 1 seul nbComps pour le moment => bug possible plus tard car on a sommer chaque comps
- # gradJ = <P.\partial_aF(a,ua)e_i>.
- # loop for each param
- # comp contain P, u_a is in trajState.comps[]
- #tmp_count+=1#FIXME
- for i in range(self.dS.study.nbParam):
- #read u_a
- code="self.daF.x.array[:]="+'+'.join(strdFp[i])
- exec(compile(code,'sumstring', 'exec'))
- # derivate of 2ndOrderTerm
- index = count*(2*self.dS.study.nbParam+1)+i+1
- formDpSOT= self.dS.secondOrderTerm.varfs[index]
- if (formDpSOT == None):scalarDPSOT=0
- else:
- # Assemble
- scalarDPSOT = assemble_scalar(formDpSOT)
- curF[i] = assemble_scalar(form(comp*self.daF*self.dS.dx)) + scalarDPSOT
+ # loop for each comps
+ for count,comp in enumerate(self.dS.bufAdjoint.comps):
+ #NOTE: on suppose 1 seul nbComps pour le moment => bug possible plus tard car on a sommer chaque comps
+ # gradJ = <P.\partial_aF(a,ua)e_i>.
+ # loop for each param
+ # comp contain P, u_a is in trajState.comps[]
+ #tmp_count+=1#FIXME
+ for i in range(self.dS.study.nbParam):
+ index = count*(2*self.dS.study.nbParam+1)+i+1
+ formDpSOT= self.dS.secondOrderTerm.varfs[index]
+ # test if derivate of SOT with respect of param is Null
+ if (formDpSOT == None):scalarDPSOT=0
+ else:
+ # Assemble
+ scalarDPSOT = assemble_scalar(formDpSOT)
+ # Take derivate of ST with respect of i-th param
+ formDpSourceTerm = self.varDerParam[i]
+ # test if derivate of SOT with respect of param is Null
+ if (formDpSourceTerm == None):scalarDPSourceTerm=0
+ else:
+ # Assemble
+ scalarDPSourceTerm = assemble_scalar(formDpSourceTerm)
+ curF[i] = scalarDPSourceTerm + scalarDPSOT
self.gradJ[:]+=l_dtAdj*(coef*pastF[:]+(2-coef)*curF[:])*0.5
# update of past T and past F
pastF[:]=curF[:]
diff --git a/trunk/tests/test_secndOrdrTermTxt.py b/trunk/tests/test_secndOrdrTermTxt.py
index b4cb2b81783a80c78cb861d5ff470d5bdd8dd813..02a9dcba678a7c251464f5660f3311ebd3f42216 100644
--- a/trunk/tests/test_secndOrdrTermTxt.py
+++ b/trunk/tests/test_secndOrdrTermTxt.py
@@ -1,7 +1,7 @@
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.MODELS.models import txtModel, model
from mse.system import system
from mse.simulator import simulator
from mse.timeDoOneStep import timeDoOneStep
@@ -12,17 +12,16 @@ from mse.explorTrajectory import explorTrajectory
import numpy as np
import pytest
from os import mkdir
-from os.path import isdir
+from os.path import isfile
### NOTE: without sourceTerm ? only 2nd order term read in file ?
##Â test with constant sourceTerm + other.
-
+# TODO: unit test
#un modele
nameModele="affineGrowth"
-nameSendOrdTerm="sOTConst"
sourceModele = computerEnv.modelsDir+nameModele+".txt"
nbComps=1
WITH_AFFINE_GROWTH=1
@@ -40,51 +39,62 @@ 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
+aStudy.setParam([2.,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()
+#add diffusion model txt
+
+def test_creationTxtDiffModel():
+ # aDiffModel type
+ nameSendOrdTerm="sOTConst"
+ aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
+ assert issubclass(txtDiffusionModel, model)
+
+def test_depedence():
+ nameSendOrdTerm="sOTConst"
+ aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
+ assert(ds2d1.secondOrderTerm == aDiffModel)
+ assert(aDiffModel.dS == ds2d1)
+ assert(aDiffModel.txtFile == computerEnv.modelsDir+nameSendOrdTerm)
+
+def test_valueVarfs():
+ nameSendOrdTerm="sOTConst"
+ aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
+ print("varfs = ",aDiffModel.varfs)
+ assert(aDiffModel.varfs[0]!=None)
+ assert(aDiffModel.varfs[1]!=None)
+ assert(aDiffModel.varfs[2]!=None)
+
+def test_fileCreation():
+ nameSendOrdTerm="sOTConst"
+ aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
+ assert(isfile(computerEnv.modelsDir+nameSendOrdTerm+"DVP.txt"))
+
+# Test 2
+
+def test_valueVarfs2():
+ nameSendOrdTerm="sOTConst3"
+ aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
+ print("varfs = ",aDiffModel.varfs)
+ assert(aDiffModel.varfs[0]!=None)
+ assert(aDiffModel.varfs[1]!=None)
+ assert(aDiffModel.varfs[2]==None)
+
+def test_fileCreation2():
+ nameSendOrdTerm="sOTConst3"
+ aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
+ assert(isfile(computerEnv.modelsDir+nameSendOrdTerm+"DVP.txt"))
+
+# Test 3
+
+def test_valueVarfs3():
+ nameSendOrdTerm="sOTConst4"
+ aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
+ print("varfs = ",aDiffModel.varfs)
+ assert(aDiffModel.varfs[0]!=None)
+ assert(aDiffModel.varfs[1]==None)
+ assert(aDiffModel.varfs[2]!=None)
+
+def test_fileCreation3():
+ nameSendOrdTerm="sOTConst4"
+ aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
+ assert(isfile(computerEnv.modelsDir+nameSendOrdTerm+"DVP.txt"))
diff --git a/trunk/tests/test_secndOrdrTermTxt2.py b/trunk/tests/test_secndOrdrTermTxt2.py
index de68797dc82a5b5248239d8ec473ed522045ec1e..8c0b306e724c3291fb1d089c2a75e7ab0c4c4be8 100644
--- a/trunk/tests/test_secndOrdrTermTxt2.py
+++ b/trunk/tests/test_secndOrdrTermTxt2.py
@@ -45,15 +45,10 @@ aSystem=system(aStudy)
ds2d1=dynamicalSystem(aStudy,2,"surface1.xdmf",nbComps=nbComps)
#add diffusion motion
+# param = [R, K, D]]
aStudy.setParam([1.,.5,5.])
-#aStudy.setParam([1.,.5])
aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
-#aDiffModel=diffusionModel(ds2d1)
-#aDiffModel.setDiffCoef(5.)
-#if (WITH_AFFINE_GROWTH):
m=txtModel(ds2d1,nameModele)
-# m.computeMatOnce=1
-# m.computeRhsOnce=1
aTS=None
# select the time algorithm
@@ -105,12 +100,12 @@ quaDiff=quadraticDiff(ds2d1)
# tabRef[:]=ds2d1.utm1.comps[0].x.array[:]
# exploreTrajSim.replay()
# # Compare: we will compare both value and check if they are different
-# #partial compy
+# #partial copy
# tabSim=ds2d1.utm1.comps[0].x.array
# assert(tabSim[5] != tabRef[5])
+# SIMU #
aStudy.setParam([1.2,.6,5.5])
-#aStudy.setParam([1.2,.6])
aStudy.simulator.doSimu()
### Now we will compute the gradV like in test_backward, and compare with kpp3
@@ -121,6 +116,16 @@ print("param Simu= ",aStudy.param)
quaDiff.computeGradV()
GJ = quaDiff.gradJ
print("gradV = :",GJ)
+# We comparethe result with the result of KPP3
+gradKPP3 = [0.2953702967205263, 0.11021221334867048, 0.0013629142893941196]
+normKPP3 =np.sqrt( 0.0013629142893941196**2 + 0.11021221334867048**2 + 0.2953702967205263**2 )
+normDiff = 0
+for i in range(len(gradKPP3)):
+ #NOTE we add first the smallest number => 3-i in loop
+ normDiff+= (gradKPP3[2-i] - GJ[2-i])**2
+errRelativ = np.sqrt(normDiff)
+print("errRelativ = ",errRelativ)
+assert (errRelativ < 0.05 )
#doPlot=False
#paramPath=[]