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()