From 72fd89e5160856bf045ab8e7399fbde976bdc806 Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Wed, 29 Aug 2018 17:13:37 -0400 Subject: [PATCH] MeshAdapt via Banded Interface (#788) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Resolved logging and simmetrix compilation * LOG: Turn on logging for the anisotropic mesh before and after * FUNC: started adding in a hook for true model reconstruction capability * USR: Added conditional to ensure model reconstruction only when given input flag * ERM: Cleaned up some aspects of anisotropy. * API: temporary fix for model reconstruction * Anisotropy, changed to constrained volume and compute size scale and frame in same loop * API: reconstruction function incorporated * MAINT: cleaned out some debug statements * ANISO: cleaned up some code for anisotropic size field * BLD: update hashstack for updated scorec lib * FIX: useModel incorporate into PUMIDomain * RECONSTRUCTION: Fixed bug in loading a reconstructed mesh regarding m�material arrays by adding a necessary input flag into MeshAdaptPUMIDrvr * RECONSTRUCTION: turn off old reconstruction call in NumericalSolution and set useModel variable if not in inputs * RECONSTRUCTION: enable 2D reconstruction, generalization of data structures for 2D, and add default behavior to useModel flag * RECONSTRUCTION: change the flag needed to start model reconstruction * STACK: changed default stack hash to get updated scorec pkgs * Update .hashstack_default * ANISO: added maximum aspect ratio option, defaults to 100 * MA: remove input for forced adaptation * ANISO: preserve size frame and scales post adapt * minor typo fix in doc string * save the meshes before and after adapt with step number * adapt based on solve step * Removed clamping from anisotropy, gradeMesh instead of smoothing * Changed how the adapt inputs are passed in to preserve input fields * Added feature to search for error near interface when no error target is supplied, turned clamping back on * anisotropic grading added through aspect ratio scaling. * Use the minimum height function and change the error to size relationship * error-to-size relationship fix for isotropic, setup framework for anisotropic logic * FIX: fix some errors with API * FIX: resolve test for gmshCheck by adding in targetError * API: added several inputs/parameters to MeshAdaptPUMIDrvr including hPhi, numAdaptSteps (for controlling adaptation) in NumericalSolution * FIX: avoid double counting of solve step when counting for adaptivity * API: enable transfer of epsFact into size field considerations * STATE: save state of initial implementation of parallel gradation algorithm * FEATURE: implementation of parallel gradation * FIX/API: require the use of context in order to get domain and physics * FIX: fix merge conflict * FIX: fix conflicts from merge * FIX:avoid node material conflicts by creating a field beforehand and using that as the source of classification * MAINT: modified some pList references to use context, ct // turned off goodQuality input * FIX: resolve parallel assignment of node materials through parallel communication * FIX/MAINT: cleanup of debug statements, added description, simplify nodeMaterial procedure to resolve bug * OPT: move mesh derivation routine to immediately after mesh generation * ENH: improve reliability of updateMaterial from a derived model by taking the minimum material tag * MAINT: for more consistent behavior on HPC, use Zoltan only for now * MAINT: turn off mass check since not all applications have the requisite fields * FIX: correct numbering of logging meshes to prevent overlapping names in parallel cases * MAINT: removed some old code/clutter * ENH: Added edge-based band intersection detection for interface based refinement, accounting for parallel cases as well * MAINT/DOC: cleaned up debug statements * MAINT/ENH: Added initial adapt before solve and moved some code into a function * MAINT: turn adapt back on, remove unnecessary statements after adapt in initial adapt * SAVE: state of numericalsolution with list of needed data structures to transfer/rebuild * FIX: ensure consistency in adapt workflow by setting free DOFs from RD and MCorr * FIX: matched shock capturing for VOF by transporting old dofs and performing computations * Save state for debugging, single phase vbdf time integration with adaptivity workflow is consistent * FIX: generalized postAdapt functions as member functions of relevant classes,can handle multiphase 2D vbdf case with systemStepExact=True. * Turned off vof shock capturing lagging. Need to store u^{n-2} solution to get it working properly * ENH/FIX: Added secondary initial adapt and corrected how nUpdatesTimeHistoryCalled is updated after adapt * FIX: avoid memory error by properly destorying marker field * TRY/FIX: force shock capturing lagging post adapt based on incorrect numDiff last and rectify subgrid error array for lagging * UPDATE: hard-coded parameter for interface band size * FIX: wrap initial adapt in sfconfig dependent conditional * TESTS: fix tests that require epsFact information to be passed in * TEST: didn't finish fixing test from previous commit * FIX: needed to destroy some entities to clear memory * VMS: finish parallel VMS implementation with simple all_reduce sum. Also extension of centroid computation to 3D. * FIX: Add definition of postAdaptUpdate for backward euler cfl * MAINT: add warning if users ever call postAdaptUpdate() without having a proper definition * FIX: shock capturing depends on subgrid error so sge needs to be in loop when recovering the shock capturing * FIX: resolved issues with lagging with multiple models by defining appropriate stages that represent the model stagger. Also ensured adaptive time stepping is accounted for with the dtLast modification. * MAINT: consistency in function def * FIX: MCorr no longer needs to be corrected since LS ebqe array was fixed * FIX: account for situation when dtLast does not exist * FIX: is not instead of != * FIX: populate the lagged eddy-viscosity array. This for some reason is not sufficient to yield exactly the same results between pseudoadapt and non-adaptive case. 1e-12 error. * TST: added pseudo option for adapt to enable ease of testing * TST: added gauge comparison test for post-adapt restart/solution transfer via 2D dambreak case * TST: added gauge compare test to setup.py * TEST: hack to get really crude test of gauge comparison running. will need to reformat into template of other tests * FIX/TST: gauge compare test not being run by Travis * FIX/TST: change the install path for gauge_compare test to fix Travis error * DOC: add some documentation strings to the grading functions * DOC: added some more documentation strings --- .gitignore | 1 + proteus/MeshAdaptPUMI/MeshAdaptPUMI.h | 9 +- proteus/MeshAdaptPUMI/MeshAdaptPUMI.pyx | 21 +- proteus/MeshAdaptPUMI/MeshConverter.cpp | 91 +++- proteus/MeshAdaptPUMI/MeshFields.cpp | 5 +- proteus/MeshAdaptPUMI/SizeField.cpp | 505 +++++++++++++++--- proteus/MeshAdaptPUMI/VMS.cpp | 15 +- proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp | 36 +- proteus/NumericalSolution.py | 361 +++++++++---- proteus/TimeIntegration.py | 36 ++ proteus/TransportCoefficients.py | 5 + proteus/mprans/MCorr.py | 11 + proteus/mprans/RDLS.py | 8 +- .../dambreak_Colagrossi_2D/README | 11 + .../dambreak_Colagrossi_2D/Reconstructed.dmg | 3 + .../dambreak_Colagrossi.py | 311 +++++++++++ .../dambreak_Colagrossi_so.py | 63 +++ .../dambreak_Colagrossi_2D/ls_consrv_n.py | 73 +++ .../dambreak_Colagrossi_2D/ls_consrv_p.py | 37 ++ .../dambreak_Colagrossi_2D/ls_n.py | 92 ++++ .../dambreak_Colagrossi_2D/ls_p.py | 35 ++ .../dambreak_Colagrossi_2D/redist_n.py | 94 ++++ .../dambreak_Colagrossi_2D/redist_p.py | 41 ++ .../test_MeshAdaptRestart.py | 37 ++ .../twp_navier_stokes_n.py | 107 ++++ .../twp_navier_stokes_p.py | 88 +++ .../dambreak_Colagrossi_2D/vof_n.py | 94 ++++ .../dambreak_Colagrossi_2D/vof_p.py | 45 ++ proteus/tests/MeshAdaptPUMI/test_gmshCheck.py | 6 +- .../MeshAdaptPUMI/test_poiseuilleError.py | 3 +- setup.py | 4 + 31 files changed, 1995 insertions(+), 253 deletions(-) create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/README create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/Reconstructed.dmg create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/dambreak_Colagrossi.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/dambreak_Colagrossi_so.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_consrv_n.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_consrv_p.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_n.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_p.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/redist_n.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/redist_p.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/test_MeshAdaptRestart.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/twp_navier_stokes_n.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/twp_navier_stokes_p.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/vof_n.py create mode 100644 proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/vof_p.py diff --git a/.gitignore b/.gitignore index 50a82d09d7..593a2afdf8 100644 --- a/.gitignore +++ b/.gitignore @@ -22,6 +22,7 @@ data*txt *.poly *.bak *.log.* +*.swp install_* config_* build/ diff --git a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h index ce3ec9c08a..95af4dd650 100644 --- a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h +++ b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h @@ -18,7 +18,7 @@ class MeshAdaptPUMIDrvr{ public: - MeshAdaptPUMIDrvr(double, double, int, const char*, const char*,const char*,double,double,int,double,double); + MeshAdaptPUMIDrvr(double, double, double, int, int, int, const char*, const char*,const char*,double,double,int,double,double); ~MeshAdaptPUMIDrvr(); int loadModelAndMesh(const char* modelFile, const char* meshFile); //load the model and mesh @@ -47,7 +47,7 @@ class MeshAdaptPUMIDrvr{ //Functions used to transfer information between PUMI and proteus int transferFieldToPUMI(const char* name, double const* inArray, int nVar, int nN); int transferFieldToProteus(const char* name, double* outArray, int nVar, int nN); - int transferPropertiesToPUMI(double* rho_p, double* nu_p,double* g_p, double deltaT); + int transferPropertiesToPUMI(double* rho_p, double* nu_p,double* g_p, double deltaT, double interfaceBandSize); //int transferBCtagsToProteus(int* tagArray, int idx, int* ebN, int* eN_global, double* fluxBC); //int transferBCsToProteus(); @@ -70,12 +70,15 @@ class MeshAdaptPUMIDrvr{ //Public Variables - double hmax, hmin; //bounds on mesh size + double hmax, hmin, hPhi; //bounds on mesh size int numIter; //number of iterations for MeshAdapt int nAdapt; //counter for number of adapt steps int nEstimate; //counter for number of error estimator calls int nsd; //number of spatial dimensions int maxAspect; //maximum aspect ratio + int adaptMesh; //is adaptivity turned on? + int numAdaptSteps; //Number adaptivity + double N_interface_band; //number of elements in half-band around interface double gradingFactor; //User Inputs diff --git a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.pyx b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.pyx index f452796069..9c302acf4c 100644 --- a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.pyx +++ b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.pyx @@ -17,9 +17,10 @@ cdef extern from "cmeshToolsModule.h": cdef extern from "MeshAdaptPUMI/MeshAdaptPUMI.h": cdef cppclass MeshAdaptPUMIDrvr: - MeshAdaptPUMIDrvr(double, double, int, char*, char*,char*,double,double,int,double,double) + MeshAdaptPUMIDrvr(double, double, double, int, int, int, char*, char*,char*,double,double,int,double,double) int numIter, numAdaptSteps string size_field_config, adapt_type_config + int adaptMesh int isReconstructed int loadModelAndMesh(char *, char*) int getSimmetrixBC() @@ -32,7 +33,7 @@ cdef extern from "MeshAdaptPUMI/MeshAdaptPUMI.h": int updateMaterialArrays2(Mesh&) int transferFieldToPUMI(char*, double*, int, int) int transferFieldToProteus(char*, double*, int, int) - int transferPropertiesToPUMI(double*, double*,double*,double) + int transferPropertiesToPUMI(double*, double*,double*,double,double) int transferModelInfo(int*,int*,int*,int*,int*,int*,int) int transferBCtagsToProteus(int*, int, int*, int*,double*) int transferBCsToProteus() @@ -49,13 +50,13 @@ cdef extern from "MeshAdaptPUMI/MeshAdaptPUMI.h": cdef class MeshAdaptPUMI: cdef MeshAdaptPUMIDrvr *thisptr - cdef double hmax, hmin + cdef double hmax, hmin, hPhi cdef int numIter, numAdaptSteps cdef int isReconstructed - cdef double deltaT - def __cinit__(self, hmax=100.0, hmin=1e-8, numIter=10, sfConfig="ERM",maType="test",logType="off",targetError=0,targetElementCount=0,reconstructedFlag=0,maxAspectRatio=100.0,gradingFact=1.5): + cdef int adaptMesh + def __cinit__(self, hmax=100.0, hmin=1e-8, hPhi=1e-2, adaptMesh=0, numIter=10, numAdaptSteps=10 ,sfConfig="ERM",maType="test",logType="off",targetError=0,targetElementCount=0,reconstructedFlag=0,maxAspectRatio=100.0,gradingFact=1.5): logEvent("MeshAdaptPUMI: hmax = {0} hmin = {1} numIter = {2}".format(hmax,hmin,numIter)) - self.thisptr = new MeshAdaptPUMIDrvr(hmax, hmin, numIter, sfConfig,maType,logType,targetError,targetElementCount,reconstructedFlag,maxAspectRatio,gradingFact) + self.thisptr = new MeshAdaptPUMIDrvr(hmax, hmin, hPhi, adaptMesh, numIter, numAdaptSteps, sfConfig,maType,logType,targetError,targetElementCount,reconstructedFlag,maxAspectRatio,gradingFact) def __dealloc__(self): del self.thisptr def size_field_config(self): @@ -64,6 +65,10 @@ cdef class MeshAdaptPUMI: return self.thisptr.adapt_type_config def numIter(self): return self.thisptr.numIter + def adaptMesh(self): + return self.thisptr.adaptMesh + def numAdaptSteps(self): + return self.thisptr.numAdaptSteps def isReconstructed(self): return self.thisptr.isReconstructed def loadModelAndMesh(self, geomName, meshName): @@ -99,11 +104,11 @@ cdef class MeshAdaptPUMI: def transferFieldToProteus(self, name, np.ndarray[np.double_t,ndim=2,mode="c"] outArray): outArray = np.ascontiguousarray(outArray) return self.thisptr.transferFieldToProteus(name, &outArray[0,0], outArray.shape[1], outArray.shape[0]) - def transferPropertiesToPUMI(self, np.ndarray[np.double_t,ndim=1,mode="c"] rho, np.ndarray[np.double_t,ndim=1,mode="c"] nu, np.ndarray[np.double_t,ndim=1,mode="c"] g, deltaT): + def transferPropertiesToPUMI(self, np.ndarray[np.double_t,ndim=1,mode="c"] rho, np.ndarray[np.double_t,ndim=1,mode="c"] nu, np.ndarray[np.double_t,ndim=1,mode="c"] g, double deltaT, double interfaceBandSize): rho = np.ascontiguousarray(rho) nu = np.ascontiguousarray(nu) g = np.ascontiguousarray(g) - return self.thisptr.transferPropertiesToPUMI(&rho[0],&nu[0],&g[0],deltaT) + return self.thisptr.transferPropertiesToPUMI(&rho[0],&nu[0],&g[0],deltaT,interfaceBandSize) def transferModelInfo(self, np.ndarray[int,ndim=1,mode="c"] numModelEntities, np.ndarray[int,ndim=2,mode="c"] edges, np.ndarray[int,ndim=2,mode="c"] faces, diff --git a/proteus/MeshAdaptPUMI/MeshConverter.cpp b/proteus/MeshAdaptPUMI/MeshConverter.cpp index 46473d793c..03e3409962 100644 --- a/proteus/MeshAdaptPUMI/MeshConverter.cpp +++ b/proteus/MeshAdaptPUMI/MeshConverter.cpp @@ -1,6 +1,7 @@ #include #include "MeshAdaptPUMI.h" +#include #include "mesh.h" #include @@ -513,26 +514,98 @@ int MeshAdaptPUMIDrvr::updateMaterialArrays2(Mesh& mesh) apf::MeshIterator* it; apf::MeshEntity* f; + //first associate all nodes with a material tag and synchronize fields to avoid mismatches + //The procedure is to have each vertex look for its classification. + //If it is classified in the region, then it is interior. + //Else, loop over adjacent faces and stop at first instance of mesh face classified on model boundary and take tag. + //If there are no such adjacent mesh faces, then set the value to be -1. This should only happen if the vertex is a shared entity. + //If the vertex is shared, communicate value to remote copies. + //When receiving values, if the current value is -1, write to field the received value. Otherwise, do nothing. + + apf::Field* nodeMaterials = apf::createLagrangeField(m, "nodeMaterials", apf::SCALAR, 1); + it = m->begin(0); + PCU_Comm_Begin(); + while(f = m->iterate(it)) + { + geomEnt = m->toModel(f); + //if classified in a region + if(m->getModelType(geomEnt) == m->getDimension()) + { + apf::setScalar(nodeMaterials,f,0,0); + } + else + { + apf::Adjacent vert_adjFace; + m->getAdjacent(f,m->getDimension()-1,vert_adjFace); + apf::MeshEntity* face; + for(int i =0; itoModel(face); + + //IF mesh face is classified on boundary + if(m->getModelType(geomEnt) == m->getDimension()-1) + { + geomTag = m->getModelTag(geomEnt); + apf::setScalar(nodeMaterials,f,0,geomTag); + if(m->isShared(f)) + { + apf::Copies remotes; + m->getRemotes(f,remotes); + for(apf::Copies::iterator iter = remotes.begin(); iter != remotes.end(); ++iter) + { + PCU_COMM_PACK(iter->first,iter->second); + PCU_COMM_PACK(iter->first,geomTag); + } + } + break; + } + if(i == vert_adjFace.getSize()-1 ) + apf::setScalar(nodeMaterials,f,0,-1); + } + } + } + PCU_Comm_Send(); + while(PCU_Comm_Receive()) + { + PCU_COMM_UNPACK(f); + PCU_COMM_UNPACK(geomTag); + int currentTag = apf::getScalar(nodeMaterials,f,0); + int newTag = std::min(currentTag,geomTag); + //if vertex is not interior and had no adjacent faces, take received values + //else take minimum value of all tags + if(currentTag==-1) + apf::setScalar(nodeMaterials,f,0,geomTag); + else + apf::setScalar(nodeMaterials,f,0,newTag); + } + //Ensure there are no mismatches across parts and then assign node materials + apf::synchronize(nodeMaterials); + it = m->begin(0); + while(f=m->iterate(it)) + { + int vID = localNumber(f); + mesh.nodeMaterialTypes[vID] = apf::getScalar(nodeMaterials,f,0); + } + //First iterate over all faces in 3D, get the model tag and apply to all downward adjacencies int dim = m->getDimension()-1; it = m->begin(dim); - while(f = m->iterate(it)){ + while(f = m->iterate(it)) + { int i = localNumber(f); geomEnt = m->toModel(f); geomTag = m->getModelTag(geomEnt); - if(m->getModelType(geomEnt) == dim){ + if(m->getModelType(geomEnt) == dim) + { mesh.elementBoundaryMaterialTypes[i] = geomTag; - apf::Adjacent face_adjVert; - m->getAdjacent(f,0,face_adjVert); - for(int j=0;jend(it); + apf::destroyField(nodeMaterials); std::cout<<"Finished faces and verts\n"; + //Loop over regions dim = m->getDimension(); it = m->begin(dim); @@ -555,7 +628,7 @@ int MeshAdaptPUMIDrvr::updateMaterialArrays2(Mesh& mesh) * scorec/core. This may be added into scorec/core eventually and removed. */ -#include +//#include #include "apfConvert.h" #include "apfMesh2.h" #include "apf.h" diff --git a/proteus/MeshAdaptPUMI/MeshFields.cpp b/proteus/MeshAdaptPUMI/MeshFields.cpp index fda1b16f6e..e24bd1b922 100644 --- a/proteus/MeshAdaptPUMI/MeshFields.cpp +++ b/proteus/MeshAdaptPUMI/MeshFields.cpp @@ -114,7 +114,7 @@ int MeshAdaptPUMIDrvr::transferFieldToProteus(const char* name, double* outArray return 0; } -int MeshAdaptPUMIDrvr::transferPropertiesToPUMI(double* rho_p, double* nu_p, double *g_p, double deltaT) +int MeshAdaptPUMIDrvr::transferPropertiesToPUMI(double* rho_p, double* nu_p, double *g_p, double deltaT, double interfaceBandSize) /** * @brief Transfer material properties to the MeshAdaptPUMI class * @@ -134,7 +134,8 @@ int MeshAdaptPUMIDrvr::transferPropertiesToPUMI(double* rho_p, double* nu_p, dou } else if(nsd==3){ g[0] = g_p[0]; g[1] = g_p[1]; g[2] = g_p[2]; - } + } + N_interface_band = interfaceBandSize; delta_T = deltaT; return 0; } diff --git a/proteus/MeshAdaptPUMI/SizeField.cpp b/proteus/MeshAdaptPUMI/SizeField.cpp index 4ea4c13bf1..d794831003 100644 --- a/proteus/MeshAdaptPUMI/SizeField.cpp +++ b/proteus/MeshAdaptPUMI/SizeField.cpp @@ -11,6 +11,7 @@ #include #include #include +#include static void SmoothField(apf::Field *f); void gradeAnisoMesh(apf::Mesh* m,double gradingFactor); @@ -18,92 +19,122 @@ void gradeAspectRatio(apf::Mesh* m, int idx, double gradingFactor); /* Based on the distance from the interface epsilon can be controlled to determine thickness of refinement near the interface */ -static double isotropicFormula(double phi, double dphi, double verr, double hmin, double hmax, double phi_s = 0) +static double isotropicFormula(double phi, double dphi, double verr, double hmin, double hmax, double phi_s = 0, double epsFact = 0) { double size; double dphi_size_factor; double v_size_factor; - //This is just a hack for now. This disable the refinement over phi and does it over phi_s - // if (phi_s != 0.0) - // { - if (fabs(phi_s) < 5.0 * hmin) + if (fabs(phi_s) < (epsFact*7.5) * hmin) return hmin; else return hmax; - // } - // else - // { - // if (fabs(phi) < 5.0 * hmin) - // { - // dphi_size_factor = fmax(hmin / 10.0, fmin(1.0, pow(((hmin / 1000.0) / fabs(dphi + 1.0e-8)), 1.0 / 2.0))); - // size = hmin * dphi_size_factor; - // } - // else - // size = hmax; - - // size = fmax(hmin / 100.0, fmin(size, 0.001 / (verr + 1.0e-8))); - - // return size; - // } } +static void setSizeField(apf::Mesh2 *m,apf::MeshEntity *vertex,double h,apf::MeshTag *marker,apf::Field* sizeField) +//helper function for banded interface to facilitate with setting the proper mesh size and parallel communication +{ + int isMarked=0; + if(m->hasTag(vertex,marker)) + isMarked=1; + double h_new; + if(isMarked) + h_new = std::min(h,apf::getScalar(sizeField,vertex,0)); + else + { + h_new = h; + int newMark = 1; + m->setIntTag(vertex,marker,&newMark); + } + apf::setScalar(sizeField,vertex,0,h_new); + + //Parallel Communication with owning copy + if(!m->isOwned(vertex)) + { + apf::Copies remotes; + m->getRemotes(vertex,remotes); + int owningPart=m->getOwner(vertex); + PCU_COMM_PACK(owningPart, remotes[owningPart]); + PCU_COMM_PACK(owningPart, h_new); + } +} + + int MeshAdaptPUMIDrvr::calculateSizeField() +//Implementation of banded interface, edge intersection algorithm +//If mesh edge intersects the 0 level-set, then the adjacent edges need to be refined { freeField(size_iso); size_iso = apf::createLagrangeField(m, "proteus_size", apf::SCALAR, 1); - apf::MeshIterator *it = m->begin(0); - apf::MeshEntity *v; apf::Field *phif = m->findField("phi"); assert(phif); - //////////////////////////////////////// - apf::Field *phisError = m->findField("phi_s"); - assert(phisError); - ///////////////////////////////////////// - apf::Field *phiCorr = m->findField("phiCorr"); - assert(phiCorr); - apf::Field *velocityError = m->findField("velocityError"); - assert(phiCorr); - while ((v = m->iterate(it))) + + apf::MeshTag* vertexMarker = m->createIntTag("vertexMarker",1); + apf::MeshIterator *it = m->begin(1); + apf::MeshEntity *edge; + + double safetyFactor = 2.0; //need to make this user defined + double L_band = N_interface_band*hPhi*safetyFactor; + + PCU_Comm_Begin(); + while ((edge = m->iterate(it))) { - double phi = apf::getScalar(phif, v, 0); - double phi_s = apf::getScalar(phisError, v, 0); - // double dphi = apf::getScalar(phiCorr, v, 0); - // double verr = apf::getScalar(velocityError, v, 0); - double size = isotropicFormula(0.0, 0.0, 0.0, hmin, hmax, phi_s); - apf::setScalar(size_iso, v, 0, size); - } - m->end(it); - /* - If you just smooth then hmax will just diffuse into the hmin band - and you won't really get a band around phi=0 with uniform diameter - hmin. Instead, reset to hmin after each smooth within the band in - order to ensure the band uses hmin. Iterate on that process until - changes in the smoothed size are less than 50% of hmin. - */ -/* - double err_h_max=hmax; - //int its=0; - //while (err_h_max > 0.5*hmin && its < 200) - for(int its=0;its<200;its++) + apf::Adjacent edge_adjVerts; + m->getAdjacent(edge,0,edge_adjVerts); + apf::MeshEntity *vertex1 = edge_adjVerts[0]; + apf::MeshEntity *vertex2 = edge_adjVerts[1]; + double phi1 = apf::getScalar(phif,vertex1,0); + double phi2 = apf::getScalar(phif,vertex2,0); + int caseNumber = 1; + if(std::fabs(phi1)>L_band) + caseNumber++; + if(std::fabs(phi2)>L_band) + caseNumber++; + + if(caseNumber==1 || caseNumber == 2) + { + setSizeField(m,vertex1,hPhi,vertexMarker,size_iso); + setSizeField(m,vertex2,hPhi,vertexMarker,size_iso); + } + else { - its++; - SmoothField(size_iso); - err_h_max=0.0; - it = m->begin(0); - while ((v = m->iterate(it))) { - double phi = apf::getScalar(phif, v, 0); - double dphi = apf::getScalar(phiCorr, v, 0); - double verr = apf::getScalar(velocityError, v, 0); - double size_current = apf::getScalar(size_iso, v, 0); - double size = fmin(size_current,isotropicFormula(phi, dphi, verr, hmin, hmax)); - err_h_max = fmax(err_h_max,fabs(size_current-size)); - apf::setScalar(size_iso, v, 0, size); + if (phi1*phi2 <0) + { + setSizeField(m,vertex1,hPhi,vertexMarker,size_iso); + setSizeField(m,vertex2,hPhi,vertexMarker,size_iso); + } + else + { + setSizeField(m,vertex1,hmax,vertexMarker,size_iso); + setSizeField(m,vertex2,hmax,vertexMarker,size_iso); } - m->end(it); } - PCU_Barrier(); -*/ + + }//end while + + PCU_Comm_Send(); + + //Take minimum between received value and current value + apf::MeshEntity *ent; + double h_received; + while(PCU_Comm_Receive()) + { + //Note: the only receiving entities should be owning copies + PCU_COMM_UNPACK(ent); + PCU_COMM_UNPACK(h_received); + //take minimum of received values + double h_current = apf::getScalar(size_iso,ent,0); + double h_final = std::min(h_current,h_received); + apf::setScalar(size_iso,ent,0,h_final); + } + + //Synchronization has all remote copies track the owning copy value + apf::synchronize(size_iso); + m->end(it); + + //Grade the Mesh gradeMesh(); + + m->destroyTag(vertexMarker); return 0; } @@ -719,6 +750,7 @@ void getTargetError(apf::Mesh* m, apf::Field* errField, double &target_error,dou element = apf::createMeshElement(m, ent); vofElem = apf::createElement(interfaceField,element); double vofVal = apf::getScalar(vofElem,apf::Vector3(1./3.,1./3.,1./3.)); + if(vofVal < 0.9 && vofVal > 0.1){ //at the interface double errorValue = apf::getScalar(errField,ent,0); errVect.push_back(errorValue); @@ -799,7 +831,6 @@ int MeshAdaptPUMIDrvr::getERMSizeField(double err_total) target_error = err_total/sqrt(m->count(nsd)); } - // Get domain volume // should only need to be computed once unless geometry is complex if (domainVolume == 0) @@ -1033,7 +1064,7 @@ int MeshAdaptPUMIDrvr::testIsotropicSizeField() return 0; } -void gradeSizeModify(apf::Mesh* m, double gradingFactor, +int gradeSizeModify(apf::Mesh* m, double gradingFactor, double size[2], apf::Adjacent edgAdjVert, apf::Adjacent vertAdjEdg, std::queue &markedEdges, @@ -1041,6 +1072,8 @@ void gradeSizeModify(apf::Mesh* m, double gradingFactor, int fieldType, int vecPos, //which idx of sizeVec to modify int idxFlag) + +//General function to actually modify sizes { //Determine a switching scheme depending on which vertex needs a modification int idx1,idx2; @@ -1055,22 +1088,39 @@ void gradeSizeModify(apf::Mesh* m, double gradingFactor, int marker[3] = {0,1,0}; double marginVal = 0.01; + int needsParallel=0; if(fieldType == apf::SCALAR){ apf::Field* size_iso = m->findField("proteus_size"); - if(size[idx1]>gradingFactor*size[idx2]){ - size[idx1] = gradingFactor*size[idx2]; - apf::setScalar(size_iso,edgAdjVert[idx1],0,size[idx1]); - m->getAdjacent(edgAdjVert[idx1], 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - //if edge is not already marked - if(!marker[2]){ - m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]); - markedEdges.push(vertAdjEdg[i]); + + if(size[idx1]>(gradingFactor*size[idx2])*(1+marginVal)) + { + if(m->isOwned(edgAdjVert[idx1])) + { + size[idx1] = gradingFactor*size[idx2]; + apf::setScalar(size_iso,edgAdjVert[idx1],0,size[idx1]); + m->getAdjacent(edgAdjVert[idx1], 1, vertAdjEdg); + for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); + //if edge is not already marked + if(!marker[2]){ + m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]); + markedEdges.push(vertAdjEdg[i]); + } } + } //end isOwned + else + { //Pack information to owning processor + needsParallel=1; + apf::Copies remotes; + m->getRemotes(edgAdjVert[idx1],remotes); + double newSize = gradingFactor*size[idx2]; + int owningPart=m->getOwner(edgAdjVert[idx1]); + PCU_COMM_PACK(owningPart, remotes[owningPart]); + PCU_COMM_PACK(owningPart,newSize); } } + }//end if apf::SCALAR else{ apf::Field* size_scale = m->findField("proteus_size_scale"); @@ -1096,23 +1146,227 @@ void gradeSizeModify(apf::Mesh* m, double gradingFactor, } } } + return needsParallel; +} + +void markEdgesInitial(apf::Mesh* m, std::queue &markedEdges,double gradingFactor) +//Function used to initially determine which edges need to be considered for gradation +{ + //marker structure for 0) not marked 1) marked 2)storage + int marker[3] = {0,1,0}; + + double size[2]; + apf::MeshTag* isMarked = m->findTag("isMarked"); + apf::Field* size_iso = m->findField("proteus_size"); + apf::Adjacent edgAdjVert; + apf::MeshEntity* edge; + apf::MeshIterator* it = m->begin(1); + while((edge=m->iterate(it))){ + m->getAdjacent(edge, 0, edgAdjVert); + for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ + size[i]=apf::getScalar(size_iso,edgAdjVert[i],0); + } + if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){ + //add edge to a queue + markedEdges.push(edge); + //tag edge to indicate that it is part of queue + m->setIntTag(edge,isMarked,&marker[1]); + } + else{ + m->setIntTag(edge,isMarked,&marker[0]); + } + } + m->end(it); +} + +int serialGradation(apf::Mesh* m, std::queue &markedEdges,double gradingFactor) +//Function used loop over the mesh edge queue for gradation and modify the sizes +{ + double size[2]; + //marker structure for 0) not marked 1) marked 2)storage + int marker[3] = {0,1,0}; + apf::MeshTag* isMarked = m->findTag("isMarked"); + apf::Field* size_iso = m->findField("proteus_size"); + apf::Adjacent edgAdjVert; + apf::Adjacent vertAdjEdg; + apf::MeshEntity* edge; + apf::MeshIterator* it = m->begin(1); + int needsParallel=0; + + //perform serial gradation while packing necessary info for parallel + while(!markedEdges.empty()){ + edge = markedEdges.front(); + m->getAdjacent(edge, 0, edgAdjVert); + for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ + size[i] = apf::getScalar(size_iso,edgAdjVert[i],0); + } + + needsParallel+=gradeSizeModify(m, gradingFactor, size, edgAdjVert, + vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 0); + needsParallel+=gradeSizeModify(m, gradingFactor, size, edgAdjVert, + vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 1); + + m->setIntTag(edge,isMarked,&marker[0]); + markedEdges.pop(); + } + return needsParallel; } int MeshAdaptPUMIDrvr::gradeMesh() //Function to grade isotropic mesh through comparison of edge vertex size ratios -//For simplicity, we do not bother with accounting for entities across partitions +//This implementation accounts for parallel meshes as well +//First do serial gradation. +//If a shared entity has its size modified, then send new size to owning copy. +//After full loop over entities, have owning copy take minimum of all sizes received +//Flag adjacent entities to owning copy. +//Communicate to remote copies that a size was modified, and flag adjacent edges to remote copies for further gradation { // if(comm_rank==0) std::cout<<"Starting grading\n"; + apf::MeshEntity* edge; + apf::Adjacent edgAdjVert; + apf::Adjacent vertAdjEdg; + double size[2]; + std::queue markedEdges; + apf::MeshTag* isMarked = m->createIntTag("isMarked",1); + + //marker structure for 0) not marked 1) marked 2)storage + int marker[3] = {0,1,0}; + + apf::MeshIterator* it; + markEdgesInitial(m,markedEdges,gradingFactor); + + int needsParallel=1; + int nCount=1; + while(needsParallel) + { + PCU_Comm_Begin(); + needsParallel = serialGradation(m,markedEdges,gradingFactor); + + PCU_Add_Ints(&needsParallel,1); + if(comm_rank==0) + std::cerr<<"Sending size info for gradation"< updateRemoteVertices; + + apf::Copies remotes; + //owning copies are receiving + while(PCU_Comm_Receive()) + { + PCU_COMM_UNPACK(ent); + PCU_COMM_UNPACK(receivedSize); + + if(!m->isOwned(ent)){ + std::cout<<"THERE WAS AN ERROR"<getAdjacent(ent, 1, vertAdjEdg); + for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); + if(!marker[2]) + { + markedEdges.push(edge); + //tag edge to indicate that it is part of queue + m->setIntTag(edge,isMarked,&marker[1]); + } + } + updateRemoteVertices.push(ent); + } + + PCU_Comm_Begin(); + + while(!updateRemoteVertices.empty()) + { + ent = updateRemoteVertices.front(); + //get remote copies and send updated mesh sizes + m->getRemotes(ent,remotes); + currentSize = apf::getScalar(size_iso,ent,0); + for(apf::Copies::iterator iter=remotes.begin(); iter!=remotes.end();++iter) + { + PCU_COMM_PACK(iter->first, iter->second); + } + updateRemoteVertices.pop(); + } + + PCU_Comm_Send(); + //while remote copies are receiving + while(PCU_Comm_Receive()) + { + //unpack + PCU_COMM_UNPACK(ent); + //PCU_COMM_UNPACK(receivedSize); + assert(!m->isOwned(ent)); + + if(m->isOwned(ent)){ + std::cout<<"Problem occurred\n"; + std::exit(1); + } + + //add adjacent edges into Q + m->getAdjacent(ent, 1, vertAdjEdg); + for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); + if(!marker[2]) + { + markedEdges.push(edge); + //tag edge to indicate that it is part of queue + m->setIntTag(edge,isMarked,&marker[1]); + } + } + } + apf::synchronize(size_iso); + + } //end outer while + + //Cleanup of edge marker field + it = m->begin(1); + while((edge=m->iterate(it))){ + m->removeTag(edge,isMarked); + } + m->end(it); + m->destroyTag(isMarked); + + //apf::synchronize(size_iso); + if(comm_rank==0) + std::cout<<"Completed grading\n"; + return needsParallel; +} + +void gradeAnisoMesh(apf::Mesh* m) +//Function to grade anisotropic mesh through comparison of edge vertex aspect ratios and minimum sizes +//For simplicity, we do not bother with accounting for entities across partitions +{ + // + //if(comm_rank==0) + // std::cout<<"Starting anisotropic grading\n"; apf::MeshIterator* it = m->begin(1); apf::MeshEntity* edge; apf::Adjacent edgAdjVert; apf::Adjacent vertAdjEdg; - //double gradingFactor = 1.5; + double gradingFactor = 1.3; double size[2]; + apf::Vector3 sizeVec; std::queue markedEdges; apf::MeshTag* isMarked = m->createIntTag("isMarked",1); + apf::Field* size_scale = m->findField("proteus_size_scale"); //marker structure for 0) not marked 1) marked 2)storage int marker[3] = {0,1,0}; @@ -1120,13 +1374,15 @@ int MeshAdaptPUMIDrvr::gradeMesh() while((edge=m->iterate(it))){ m->getAdjacent(edge, 0, edgAdjVert); for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ - size[i]=apf::getScalar(size_iso,edgAdjVert[i],0); + apf::getVector(size_scale,edgAdjVert[i],0,sizeVec); + size[i]=sizeVec[0]; } if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){ //add edge to a queue markedEdges.push(edge); //tag edge to indicate that it is part of queue m->setIntTag(edge,isMarked,&marker[1]); + } else{ m->setIntTag(edge,isMarked,&marker[0]); @@ -1137,12 +1393,20 @@ int MeshAdaptPUMIDrvr::gradeMesh() edge = markedEdges.front(); m->getAdjacent(edge, 0, edgAdjVert); for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ - size[i] = apf::getScalar(size_iso,edgAdjVert[i],0); + apf::getVector(size_scale,edgAdjVert[i],0,sizeVec); + size[i]=sizeVec[0]; } - //This can be simplified with one function call + gradeSizeModify(m, gradingFactor, size, edgAdjVert, + vertAdjEdg, markedEdges, isMarked, apf::VECTOR,0, 0); + gradeSizeModify(m, gradingFactor, size, edgAdjVert, + vertAdjEdg, markedEdges, isMarked, apf::VECTOR,0, 1); + +/* if(size[0]>gradingFactor*size[1]){ size[0] = gradingFactor*size[1]; - apf::setScalar(size_iso,edgAdjVert[0],0,size[0]); + apf::getVector(size_scale,edgAdjVert[0],0,sizeVec); + sizeVec[0] = size[0]; + apf::setVector(size_scale,edgAdjVert[0],0,sizeVec); m->getAdjacent(edgAdjVert[0], 1, vertAdjEdg); for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); @@ -1155,7 +1419,9 @@ int MeshAdaptPUMIDrvr::gradeMesh() } if(size[1]>gradingFactor*size[0]){ size[1] = gradingFactor*size[0]; - apf::setScalar(size_iso,edgAdjVert[1],0,size[1]); + apf::getVector(size_scale,edgAdjVert[1],0,sizeVec); + sizeVec[0] = size[1]; + apf::setVector(size_scale,edgAdjVert[1],0,sizeVec); m->getAdjacent(edgAdjVert[1], 1, vertAdjEdg); for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); @@ -1166,19 +1432,82 @@ int MeshAdaptPUMIDrvr::gradeMesh() } } } +*/ m->setIntTag(edge,isMarked,&marker[0]); markedEdges.pop(); } + it = m->begin(1); + while((edge=m->iterate(it))){ + m->removeTag(edge,isMarked); + } + m->end(it); + m->destroyTag(isMarked); + apf::synchronize(size_scale); + //if(comm_rank==0) + // std::cout<<"Completed minimum size grading\n"; +} +void gradeAspectRatio(apf::Mesh* m,int idx) +//Function to grade anisotropic mesh through comparison of edge vertex aspect ratios and minimum sizes +//For simplicity, we do not bother with accounting for entities across partitions +{ + std::cout<<"Entered function\n"; + apf::MeshIterator* it = m->begin(1); + apf::MeshEntity* edge; + apf::Adjacent edgAdjVert; + apf::Adjacent vertAdjEdg; + double gradingFactor = 1.3; + double size[2]; + apf::Vector3 sizeVec; + std::queue markedEdges; + apf::MeshTag* isMarked = m->createIntTag("isMarked",1); + apf::Field* size_scale = m->findField("proteus_size_scale"); + + //marker structure for 0) not marked 1) marked 2)storage + int marker[3] = {0,1,0}; + + while((edge=m->iterate(it))){ + m->getAdjacent(edge, 0, edgAdjVert); + for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ + apf::getVector(size_scale,edgAdjVert[i],0,sizeVec); + size[i]=sizeVec[idx]/sizeVec[0]; + } + if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){ + //add edge to a queue + markedEdges.push(edge); + //tag edge to indicate that it is part of queue + m->setIntTag(edge,isMarked,&marker[1]); + + } + else{ + m->setIntTag(edge,isMarked,&marker[0]); + } + } + m->end(it); + + std::cout<<"Got queue of size "<getAdjacent(edge, 0, edgAdjVert); + for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ + apf::getVector(size_scale,edgAdjVert[i],0,sizeVec); + size[i]=sizeVec[idx]/sizeVec[0]; + } + gradeSizeModify(m, gradingFactor, size, edgAdjVert, + vertAdjEdg, markedEdges, isMarked, apf::VECTOR, idx, 0); + gradeSizeModify(m, gradingFactor, size, edgAdjVert, + vertAdjEdg, markedEdges, isMarked, apf::VECTOR, idx, 1); + + m->setIntTag(edge,isMarked,&marker[0]); + markedEdges.pop(); + } it = m->begin(1); while((edge=m->iterate(it))){ m->removeTag(edge,isMarked); } m->end(it); m->destroyTag(isMarked); - apf::synchronize(size_iso); - if(comm_rank==0) - std::cout<<"Completed grading\n"; + apf::synchronize(size_scale); } void gradeAnisoMesh(apf::Mesh* m,double gradingFactor) diff --git a/proteus/MeshAdaptPUMI/VMS.cpp b/proteus/MeshAdaptPUMI/VMS.cpp index a674a41023..8c925c3c1b 100644 --- a/proteus/MeshAdaptPUMI/VMS.cpp +++ b/proteus/MeshAdaptPUMI/VMS.cpp @@ -248,9 +248,9 @@ void MeshAdaptPUMIDrvr::get_VMS_error(double &total_error) */ //H1 error compute nu_err at centroid and compute residual - qpt[0] = 1./3.; - qpt[1] = 1./3.; - qpt[2] = 0.0; + qpt[0] = 0.0; qpt[1] = 0.0; qpt[2] = 0.0; //ensure it is initialized to 0 + for(int i=0; i0) - std::cout<<"Improper parallel implementation\n"; + PCU_Add_Doubles(&VMSerrTotalL2,1); + PCU_Add_Doubles(&VMSerrTotalH1,1); if(PCU_Comm_Self()==0) std::cout<writeNative(namebuffer); - /* Code to output size scale and frame apf::MeshIterator* it = m->begin(0); apf::MeshEntity* test; @@ -453,14 +453,15 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh() } ma::validateInput(in); in->shouldRunPreZoltan = true; - in->shouldRunMidParma = true; - in->shouldRunPostParma = true; + in->shouldRunMidZoltan = true; + in->shouldRunPostZoltan = true; + //in->shouldRunMidParma = true; + //in->shouldRunPostParma = true; in->maximumImbalance = 1.05; in->maximumIterations = numIter; in->shouldSnap = false; - in->goodQuality = 0.008; - - double mass_before = getTotalMass(); + //in->goodQuality = 0.16;//0.027; + //double mass_before = getTotalMass(); double t1 = PCU_Time(); //ma::adapt(in); @@ -468,9 +469,9 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh() double t2 = PCU_Time(); m->verify(); - double mass_after = getTotalMass(); - PCU_Add_Doubles(&mass_before,1); - PCU_Add_Doubles(&mass_after,1); + //double mass_after = getTotalMass(); + //PCU_Add_Doubles(&mass_before,1); + //PCU_Add_Doubles(&mass_after,1); if(comm_rank==0 && logging_config=="on"){ /* std::ios::fmtflags saved(std::cout.flags()); @@ -481,21 +482,22 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh() myfile.open("adapt_timing.txt", std::ios::app); myfile << t2-t1<writeNative(namebuffer); } //isReconstructed = 0; //this is needed to maintain consistency with the post-adapt conversion back to Proteus diff --git a/proteus/NumericalSolution.py b/proteus/NumericalSolution.py index 1ec3eb2408..a73063171e 100644 --- a/proteus/NumericalSolution.py +++ b/proteus/NumericalSolution.py @@ -508,34 +508,6 @@ def __init__(self,so,pList,nList,sList,opts,simFlagsList=None): nLayersOfOverlap=n.nLayersOfOverlapForParallel, parallelPartitioningType=n.parallelPartitioningType) - if hasattr(p.domain,"PUMIMesh") and not isinstance(p.domain,Domain.PUMIDomain) : - logEvent("Reconstruct based on Proteus, convert PUMI mesh to Proteus") - p = self.pList[0] - n = self.nList[0] - - from scipy import spatial - meshVertexTree = spatial.cKDTree(mesh.nodeArray) - meshVertex2Model= [0]*mesh.nNodes_owned - for idx,vertex in enumerate(self.pList[0].domain.vertices): - if(self.pList[0].nd==2 and len(vertex) == 2): #there might be a smarter way to do this - vertex.append(0.0) #need to make a 3D coordinate - closestVertex = meshVertexTree.query(vertex) - meshVertex2Model[closestVertex[1]] = 1 - - isModelVert = numpy.asarray(meshVertex2Model).astype("i") - - meshBoundaryConnectivity = numpy.zeros((mesh.nExteriorElementBoundaries_global,2+pList[0].nd),dtype=numpy.int32) - for elementBdyIdx in range(len(mesh.exteriorElementBoundariesArray)): - exteriorIdx = mesh.exteriorElementBoundariesArray[elementBdyIdx] - meshBoundaryConnectivity[elementBdyIdx][0] = mesh.elementBoundaryMaterialTypes[exteriorIdx] - meshBoundaryConnectivity[elementBdyIdx][1] = mesh.elementBoundaryElementsArray[exteriorIdx][0] - meshBoundaryConnectivity[elementBdyIdx][2] = mesh.elementBoundaryNodesArray[exteriorIdx][0] - meshBoundaryConnectivity[elementBdyIdx][3] = mesh.elementBoundaryNodesArray[exteriorIdx][1] - if(self.pList[0].nd==3): - meshBoundaryConnectivity[elementBdyIdx][4] = mesh.elementBoundaryNodesArray[exteriorIdx][2] - - p.domain.PUMIMesh.reconstructFromProteus2(mesh.cmesh,isModelVert,meshBoundaryConnectivity) - mlMesh_nList.append(mlMesh) if opts.viewMesh: logEvent("Attempting to visualize mesh") @@ -560,6 +532,38 @@ def __init__(self,so,pList,nList,sList,opts,simFlagsList=None): viewBoundaryMaterialTypes=True) except: logEvent("NumericalSolution ViewMesh failed for mesh level %s" % l) + + theMesh = mlMesh.meshList[0].subdomainMesh + pCT = self.pList[0]#self.pList[0].ct + nCT = self.nList[0]#self.nList[0].ct + theDomain = pCT.domain + + if hasattr(theDomain,"PUMIMesh") and not isinstance(theDomain,Domain.PUMIDomain) : + logEvent("Reconstruct based on Proteus, convert PUMI mesh to Proteus") + + from scipy import spatial + meshVertexTree = spatial.cKDTree(theMesh.nodeArray) + meshVertex2Model= [0]*theMesh.nNodes_owned + for idx,vertex in enumerate(theDomain.vertices): + if(pCT.nd==2 and len(vertex) == 2): #there might be a smarter way to do this + vertex.append(0.0) #need to make a 3D coordinate + closestVertex = meshVertexTree.query(vertex) + meshVertex2Model[closestVertex[1]] = 1 + + isModelVert = numpy.asarray(meshVertex2Model).astype("i") + + meshBoundaryConnectivity = numpy.zeros((theMesh.nExteriorElementBoundaries_global,2+pCT.nd),dtype=numpy.int32) + for elementBdyIdx in range(len(theMesh.exteriorElementBoundariesArray)): + exteriorIdx = theMesh.exteriorElementBoundariesArray[elementBdyIdx] + meshBoundaryConnectivity[elementBdyIdx][0] = theMesh.elementBoundaryMaterialTypes[exteriorIdx] + meshBoundaryConnectivity[elementBdyIdx][1] = theMesh.elementBoundaryElementsArray[exteriorIdx][0] + meshBoundaryConnectivity[elementBdyIdx][2] = theMesh.elementBoundaryNodesArray[exteriorIdx][0] + meshBoundaryConnectivity[elementBdyIdx][3] = theMesh.elementBoundaryNodesArray[exteriorIdx][1] + if(pCT.nd==3): + meshBoundaryConnectivity[elementBdyIdx][4] = theMesh.elementBoundaryNodesArray[exteriorIdx][2] + + pCT.domain.PUMIMesh.reconstructFromProteus2(theMesh.cmesh,isModelVert,meshBoundaryConnectivity) + if so.useOneMesh: for p in pList[1:]: mlMesh_nList.append(mlMesh) try: @@ -752,8 +756,11 @@ def allocateModels(self): Profiling.memory("MultilevelNonlinearSolver for"+p.name) def PUMI2Proteus(self,mesh): - p0 = self.pList[0] #This can probably be cleaned up somehow - n0 = self.nList[0] + #p0 = self.pList[0] #This can probably be cleaned up somehow + #n0 = self.nList[0] + p0 = self.pList[0].ct + n0 = self.nList[0].ct + logEvent("Generating %i-level mesh from PUMI mesh" % (n0.nLevels,)) if p0.domain.nd == 3: mlMesh = MeshTools.MultilevelTetrahedralMesh( @@ -862,6 +869,11 @@ def PUMI2Proteus(self,mesh): coef.vectorName, vector) for vci in range(len(coef.vectorComponents)): lm.u[coef.vectorComponents[vci]].dof[:] = vector[:,vci] + p0.domain.PUMIMesh.transferFieldToProteus( + coef.vectorName+"_old", vector) + for vci in range(len(coef.vectorComponents)): + lm.u[coef.vectorComponents[vci]].dof_last[:] = vector[:,vci] + del vector for ci in range(coef.nc): if coef.vectorComponents is None or \ @@ -870,17 +882,20 @@ def PUMI2Proteus(self,mesh): p0.domain.PUMIMesh.transferFieldToProteus( coef.variableNames[ci], scalar) lm.u[ci].dof[:] = scalar[:,0] + p0.domain.PUMIMesh.transferFieldToProteus( + coef.variableNames[ci]+"_old", scalar) + lm.u[ci].dof_last[:] = scalar[:,0] del scalar logEvent("Attaching models on new mesh to each other") for m,ptmp,mOld in zip(self.modelList, self.pList, modelListOld): for lm, lu, lr, lmOld in zip(m.levelModelList, m.uList, m.rList,mOld.levelModelList): - save_dof=[] - for ci in range(lm.coefficients.nc): - save_dof.append( lm.u[ci].dof.copy()) - lm.u[ci].dof_last = lm.u[ci].dof.copy() + #save_dof=[] + #for ci in range(lm.coefficients.nc): + # save_dof.append( lm.u[ci].dof.copy()) + # lm.u[ci].dof_last = lm.u[ci].dof.copy() lm.setFreeDOF(lu) - for ci in range(lm.coefficients.nc): - assert((save_dof[ci] == lm.u[ci].dof).all()) + #for ci in range(lm.coefficients.nc): + # assert((save_dof[ci] == lm.u[ci].dof).all()) lm.calculateSolutionAtQuadrature() lm.timeIntegration.tLast = lmOld.timeIntegration.tLast lm.timeIntegration.t = lmOld.timeIntegration.t @@ -892,11 +907,15 @@ def PUMI2Proteus(self,mesh): m.stepController.t_model = mOld.stepController.t_model m.stepController.t_model_last = mOld.stepController.t_model_last m.stepController.substeps = mOld.stepController.substeps - # logEvent("Evaluating residuals and time integration") + + #Attach models and do sample residual calculation. The results are usually irrelevant. + #What's important right now is to re-establish the relationships between data structures. + #The necessary values will be written in later. for m,ptmp,mOld in zip(self.modelList, self.pList, modelListOld): logEvent("Attaching models to model "+ptmp.name) m.attachModels(self.modelList) logEvent("Evaluating residuals and time integration") + for m,ptmp,mOld in zip(self.modelList, self.pList, modelListOld): for lm, lu, lr, lmOld in zip(m.levelModelList, m.uList, m.rList, mOld.levelModelList): lm.timeTerm=True @@ -905,12 +924,7 @@ def PUMI2Proteus(self,mesh): lm.initializeTimeHistory() lm.timeIntegration.initializeSpaceHistory() lm.getResidual(lu,lr) - # assert(lmOld.timeIntegration.tLast == lm.timeIntegration.tLast) - # assert(lmOld.timeIntegration.t == lm.timeIntegration.t) - # assert(lmOld.timeIntegration.dt == lm.timeIntegration.dt) - #lm.coefficients.evaluate(self.t_stepSequence,lm.q) - #lm.coefficients.evaluate(self.t_stepSequence,lm.ebqe) - #lm.timeIntegration.calculateElementCoefficients(lm.q) + #lm.estimate_mt() #function is empty in all models assert(m.stepController.dt_model == mOld.stepController.dt_model) assert(m.stepController.t_model == mOld.stepController.t_model) assert(m.stepController.t_model_last == mOld.stepController.t_model_last) @@ -926,14 +940,91 @@ def PUMI2Proteus(self,mesh): self.systemStepController.controllerList.append(model) self.systemStepController.maxFailures = model.stepController.maxSolverFailures self.systemStepController.choose_dt_system() - # for m,ptmp,mOld in zip(self.modelList, self.pList, modelListOld): - # for lm, lu, lr, lmOld in zip(m.levelModelList, m.uList, m.rList, mOld.levelModelList): - # assert(lmOld.timeIntegration.tLast == lm.timeIntegration.tLast) - # assert(lmOld.timeIntegration.t == lm.timeIntegration.t) - # assert(lmOld.timeIntegration.dt == lm.timeIntegration.dt) - # assert(m.stepController.dt_model == mOld.stepController.dt_model) - # assert(m.stepController.t_model == mOld.stepController.t_model) - # assert(m.stepController.t_model_last == mOld.stepController.t_model_last) + + ##This section is to correct any differences in the quadrature point field from the old model + + #Shock capturing lagging needs to be matched + + import copy + #This sections gets beta bdf right + #import pdb; pdb.set_trace() + #self.modelList[1].levelModelList[0].u_store = copy.deepcopy(self.modelList[1].levelModelList[0].u) + #self.modelList[1].levelModelList[0].u[0].dof[:] = self.modelList[1].levelModelList[0].u[0].dof_last + #self.modelList[1].levelModelList[0].calculateElementResidual() + #self.modelList[1].levelModelList[0].q[('m_last',0)][:] = self.modelList[1].levelModelList[0].q[('m_tmp',0)] + + ##this section gets numDiff right + #self.modelList[1].levelModelList[0].u[0].dof[:] = self.modelList[1].levelModelList[0].u_store[0].dof + #self.modelList[1].levelModelList[0].u[0].dof_last[:] = self.modelList[1].levelModelList[0].u_store[0].dof_last + + #self.modelList[1].levelModelList[0].calculateElementResidual() + #self.modelList[1].levelModelList[0].q[('m_last',0)][:] = self.modelList[1].levelModelList[0].q[('m_tmp',0)] + + #if(modelListOld[1].levelModelList[0].shockCapturing.nStepsToDelay is not None and modelListOld[1].levelModelList[0].shockCapturing.nSteps > modelListOld[1].levelModelList[0].shockCapturing.nStepsToDelay): + # self.modelList[1].levelModelList[0].shockCapturing.nSteps=self.modelList[1].levelModelList[0].shockCapturing.nStepsToDelay + # self.modelList[1].levelModelList[0].shockCapturing.updateShockCapturingHistory() + + + ###Details for solution transfer + #To get shock capturing lagging correct, the numDiff array needs to be computed correctly with the u^{n} solution. + #numDiff depends on the PDE residual and can depend on the subgrid error (SGE) + #the PDE residual depends on the alpha and beta_bdf terms which depend on m_tmp from u^{n-1} as well as VOF or LS fields. + #getResidual() is used to populate m_tmp, numDiff. + #The goal is therefore to populate the nodal fields with the old solution, get m_tmp properly and lagged sge properly. + #Mimic the solver stagger with a new loop to repopulate the nodal fields with u^{n} solution. This is necessary because NS relies on the u^{n-1} field for VOF/LS + + ###This loop stores the current solution (u^n) and loads in the previous timestep solution (u^{n-1}) + for m,mOld in zip(self.modelList, modelListOld): + for lm, lu, lr, lmOld in zip(m.levelModelList, m.uList, m.rList, mOld.levelModelList): + #lm.coefficients.postAdaptStep() #MCorr needs this at the moment + lm.u_store = copy.deepcopy(lm.u) + lm.dt_store = copy.deepcopy(lm.timeIntegration.dt) + for ci in range(0,lm.coefficients.nc): + lm.u[ci].dof[:] = lm.u[ci].dof_last + lm.setFreeDOF(lu) + + #All solution fields are now in state u^{n-1} + for m,mOld in zip(self.modelList, modelListOld): + for lm, lu, lr, lmOld in zip(m.levelModelList, m.uList, m.rList, mOld.levelModelList): + lm.getResidual(lu,lr) + lm.timeIntegration.postAdaptUpdate(lmOld.timeIntegration) + + if(hasattr(lm.timeIntegration,"dtLast") and lm.timeIntegration.dtLast is not None): + lm.timeIntegration.dt = lm.timeIntegration.dtLast + + #This gets the subgrid error history correct + if(modelListOld[0].levelModelList[0].stabilization.lag and modelListOld[0].levelModelList[0].stabilization.nSteps > modelListOld[0].levelModelList[0].stabilization.nStepsToDelay): + self.modelList[0].levelModelList[0].stabilization.nSteps = self.modelList[0].levelModelList[0].stabilization.nStepsToDelay + self.modelList[0].levelModelList[0].stabilization.updateSubgridErrorHistory() + + #update the eddy-viscosity history + lm.calculateAuxiliaryQuantitiesAfterStep() + + ###This loop reloads the current solution and the previous solution into proper places + for m,mOld in zip(self.modelList, modelListOld): + for lm, lu, lr, lmOld in zip(m.levelModelList, m.uList, m.rList, mOld.levelModelList): + for ci in range(0,lm.coefficients.nc): + lm.u[ci].dof[:] = lm.u_store[ci].dof + lm.u[ci].dof_last[:] = lm.u_store[ci].dof_last + lm.setFreeDOF(lu) + lm.getResidual(lu,lr) + lm.timeIntegration.postAdaptUpdate(lmOld.timeIntegration) + lm.timeIntegration.dt = lm.dt_store + + #This gets the subgrid error history correct + if(modelListOld[0].levelModelList[0].stabilization.lag and modelListOld[0].levelModelList[0].stabilization.nSteps > modelListOld[0].levelModelList[0].stabilization.nStepsToDelay): + self.modelList[0].levelModelList[0].stabilization.nSteps = self.modelList[0].levelModelList[0].stabilization.nStepsToDelay + self.modelList[0].levelModelList[0].stabilization.updateSubgridErrorHistory() + ### + + ###Shock capturing + if(lmOld.shockCapturing and lmOld.shockCapturing.nStepsToDelay is not None and lmOld.shockCapturing.nSteps > lmOld.shockCapturing.nStepsToDelay): + lm.shockCapturing.nSteps=lm.shockCapturing.nStepsToDelay + lm.shockCapturing.updateShockCapturingHistory() + + #update the eddy-viscosity history + lm.calculateAuxiliaryQuantitiesAfterStep() + if self.archiveFlag == ArchiveFlags.EVERY_SEQUENCE_STEP: #hack for archiving initial solution on adapted mesh self.tCount+=1 @@ -943,6 +1034,64 @@ def PUMI2Proteus(self,mesh): index, self.systemStepController.t_system_last+1.0e-6) + def PUMI_transferFields(self): + p0 = self.pList[0].ct + n0 = self.nList[0].ct + + logEvent("Copying coordinates to PUMI") + p0.domain.PUMIMesh.transferFieldToPUMI("coordinates", + self.modelList[0].levelModelList[0].mesh.nodeArray) + + logEvent("Copying DOF and parameters to PUMI") + for m in self.modelList: + for lm in m.levelModelList: + coef = lm.coefficients + if coef.vectorComponents is not None: + vector=numpy.zeros((lm.mesh.nNodes_global,3),'d') + for vci in range(len(coef.vectorComponents)): + vector[:,vci] = lm.u[coef.vectorComponents[vci]].dof[:] + + p0.domain.PUMIMesh.transferFieldToPUMI( + coef.vectorName, vector) + #Transfer dof_last + for vci in range(len(coef.vectorComponents)): + vector[:,vci] = lm.u[coef.vectorComponents[vci]].dof_last[:] + p0.domain.PUMIMesh.transferFieldToPUMI( + coef.vectorName+"_old", vector) + del vector + for ci in range(coef.nc): + if coef.vectorComponents is None or \ + ci not in coef.vectorComponents: + scalar=numpy.zeros((lm.mesh.nNodes_global,1),'d') + scalar[:,0] = lm.u[ci].dof[:] + p0.domain.PUMIMesh.transferFieldToPUMI( + coef.variableNames[ci], scalar) + + #Transfer dof_last + scalar[:,0] = lm.u[ci].dof_last[:] + p0.domain.PUMIMesh.transferFieldToPUMI( + coef.variableNames[ci]+"_old", scalar) + del scalar + + scalar=numpy.zeros((lm.mesh.nNodes_global,1),'d') + + del scalar + #Get Physical Parameters + #Can we do this in a problem-independent way? + rho = numpy.array([self.pList[0].ct.rho_0, + self.pList[0].ct.rho_1]) + nu = numpy.array([self.pList[0].ct.nu_0, + self.pList[0].ct.nu_1]) + g = numpy.asarray(self.pList[0].ct.g) + if(hasattr(self,"tn")): + deltaT = self.tn-self.tn_last + else: + deltaT = 0 + epsFact = p0.epsFact_density + p0.domain.PUMIMesh.transferPropertiesToPUMI(rho,nu,g,deltaT,epsFact) + del rho, nu, g, epsFact + + def PUMI_estimateError(self): """ Estimate the error using the classical element residual method by @@ -951,14 +1100,17 @@ def PUMI_estimateError(self): p0 = self.pList[0] n0 = self.nList[0] + #p0 = self.pList[0].ct + #n0 = self.nList[0].ct + adaptMeshNow = False #will need to move this to earlier when the mesh is created #from proteus.MeshAdaptPUMI import MeshAdaptPUMI - if not hasattr(p0.domain,'PUMIMesh') and not isinstance(p0.domain,Domain.PUMIDomain) and n0.adaptMesh: - import sys - if(self.comm.size()>1 and p0.domain.MeshOptions.parallelPartitioningType!=MeshTools.MeshParallelPartitioningTypes.element): - sys.exit("The mesh must be partitioned by elements and NOT nodes for adaptivity functionality. Do this with: `domain.MeshOptions.setParallelPartitioningType('element')'.") - p0.domain.PUMIMesh=n0.MeshAdaptMesh + #if not hasattr(p0.domain,'PUMIMesh') and not isinstance(p0.domain,Domain.PUMIDomain) and p0.domain.PUMIMesh.adaptMesh(): + # import sys + # if(self.comm.size()>1 and p0.domain.MeshOptions.parallelPartitioningType!=MeshTools.MeshParallelPartitioningTypes.element): + # sys.exit("The mesh must be partitioned by elements and NOT nodes for adaptivity functionality. Do this with: `domain.MeshOptions.setParallelPartitioningType('element')'.") + # p0.domain.PUMIMesh=n0.MeshAdaptMesh #p0.domain.hasModel = n0.useModel #numModelEntities = numpy.array([len(p0.domain.vertices),len(p0.domain.segments),len(p0.domain.facets),len(p0.domain.regions)]).astype("i") ##force initialization of arrays to enable passage through to C++ code @@ -1027,12 +1179,9 @@ def PUMI_estimateError(self): #p0.domain.PUMIMesh.transferModelInfo(numModelEntities,segmentList,newFacetList,mesh2Model_v,mesh2Model_e,mesh2Model_b) #p0.domain.PUMIMesh.reconstructFromProteus(self.modelList[0].levelModelList[0].mesh.cmesh,self.modelList[0].levelModelList[0].mesh.globalMesh.cmesh,p0.domain.hasModel) if (hasattr(p0.domain, 'PUMIMesh') and - n0.adaptMesh and + p0.domain.PUMIMesh.adaptMesh() and self.so.useOneMesh and - self.nSolveSteps%n0.adaptMesh_nSteps==0): - logEvent("Copying coordinates to PUMI") - p0.domain.PUMIMesh.transferFieldToPUMI("coordinates", - self.modelList[0].levelModelList[0].mesh.nodeArray) + self.nSolveSteps%p0.domain.PUMIMesh.numAdaptSteps()==0): if (p0.domain.PUMIMesh.size_field_config() == "isotropicProteus"): p0.domain.PUMIMesh.transferFieldToPUMI("proteus_size", self.modelList[0].levelModelList[0].mesh.size_field) @@ -1051,49 +1200,8 @@ def PUMI_estimateError(self): self.modelList[0].levelModelList[0].mesh.size_scale p0.domain.PUMIMesh.transferFieldToPUMI("proteus_sizeScale", self.modelList[0].levelModelList[0].mesh.size_scale) p0.domain.PUMIMesh.transferFieldToPUMI("proteus_sizeFrame", self.modelList[0].levelModelList[0].mesh.size_frame) - p0.domain.PUMIMesh.transferFieldToPUMI("coordinates", - self.modelList[0].levelModelList[0].mesh.nodeArray) - logEvent("Copying DOF and parameters to PUMI") - for m in self.modelList: - for lm in m.levelModelList: - coef = lm.coefficients - if coef.vectorComponents is not None: - vector=numpy.zeros((lm.mesh.nNodes_global,3),'d') - for vci in range(len(coef.vectorComponents)): - vector[:,vci] = lm.u[coef.vectorComponents[vci]].dof[:] - p0.domain.PUMIMesh.transferFieldToPUMI( - coef.vectorName, vector) - #Transfer dof_last - for vci in range(len(coef.vectorComponents)): - vector[:,vci] = lm.u[coef.vectorComponents[vci]].dof_last[:] - p0.domain.PUMIMesh.transferFieldToPUMI( - "velocity_old", vector) - del vector - for ci in range(coef.nc): - if coef.vectorComponents is None or \ - ci not in coef.vectorComponents: - scalar=numpy.zeros((lm.mesh.nNodes_global,1),'d') - scalar[:,0] = lm.u[ci].dof[:] - p0.domain.PUMIMesh.transferFieldToPUMI( - coef.variableNames[ci], scalar) - del scalar - - scalar=numpy.zeros((lm.mesh.nNodes_global,1),'d') - - del scalar - #Get Physical Parameters - #Can we do this in a problem-independent way? - rho = numpy.array([self.pList[0].rho_0, - self.pList[0].rho_1]) - nu = numpy.array([self.pList[0].nu_0, - self.pList[0].nu_1]) - g = numpy.asarray(self.pList[0].g) - - deltaT_test = self.systemStepController.dt_system - deltaT = self.systemStepController.t_system-self.systemStepController.t_system_last - print("deltaT %f",deltaT," dt_system ", deltaT_test) - p0.domain.PUMIMesh.transferPropertiesToPUMI(rho,nu,g,deltaT) - del rho, nu, g, deltaT + + self.PUMI_transferFields() logEvent("Estimate Error") sfConfig = p0.domain.PUMIMesh.size_field_config() @@ -1147,11 +1255,15 @@ def PUMI_adaptMesh(self): # # self.modelList[0].levelModelList[0].numericalFlux.mesh.exteriorElementBoundariesArray, # # self.modelList[0].levelModelList[0].numericalFlux.mesh.elementBoundaryElementsArray, # # diff_flux) - p0 = self.pList[0] - n0 = self.nList[0] + + p0 = self.pList[0].ct + n0 = self.nList[0].ct sfConfig = p0.domain.PUMIMesh.size_field_config() logEvent("h-adapt mesh by calling AdaptPUMIMesh") - p0.domain.PUMIMesh.adaptPUMIMesh() + if(sfConfig=="pseudo"): + print "Testing solution transfer and restart feature of adaptation. No actual mesh adaptation!" + else: + p0.domain.PUMIMesh.adaptPUMIMesh() #code to suggest adapting until error is reduced; #not fully baked and can lead to infinite loops of adaptation @@ -1176,6 +1288,7 @@ def PUMI_adaptMesh(self): p0.domain.regList, parallel = self.comm.size() > 1, dim = p0.domain.nd) + self.PUMI2Proteus(mesh) ##chitak end Adapt @@ -1287,6 +1400,8 @@ def calculateSolution(self,runName): logEvent("Hotstarting, new tnList is"+`self.tnList`) else: self.tCount=0#time step counter + + logEvent("Attaching models and running spin-up step if requested") self.firstStep = True ##\todo get rid of firstStep flag in NumericalSolution if possible? spinup = [] @@ -1337,8 +1452,10 @@ def calculateSolution(self,runName): logEvent("Spin-up step exact called for model %s" % (m.name,),level=3) m.stepController.stepExact_model(self.tnList[1]) logEvent("Spin-Up Initializing time history for model step controller") + m.stepController.initializeTimeHistory() m.stepController.setInitialGuess(m.uList,m.rList) + solverFailed = m.solver.solveMultilevel(uList=m.uList, rList=m.rList, par_uList=m.par_uList, @@ -1407,6 +1524,28 @@ def calculateSolution(self,runName): systemStepFailed=False stepFailed=False + #### Perform an initial adapt after applying initial conditions #### + # The initial adapt is based on interface, but will eventually be generalized to any sort of initialization + # Needs to be placed here at this time because of the post-adapt routine requirements + + if (hasattr(self.pList[0].domain, 'PUMIMesh') and + self.pList[0].domain.PUMIMesh.adaptMesh() and + (self.pList[0].domain.PUMIMesh.size_field_config() == "combined" or self.pList[0].domain.PUMIMesh.size_field_config() == "interface") and + self.so.useOneMesh): + + self.PUMI_transferFields() + logEvent("Initial Adapt before Solve") + self.PUMI_adaptMesh() + + for index,p,n,m,simOutput in zip(range(len(self.modelList)),self.pList,self.nList,self.modelList,self.simOutputList): + if p.initialConditions is not None: + logEvent("Setting initial conditions for "+p.name) + m.setInitialConditions(p.initialConditions,self.tnList[0]) + + self.PUMI_transferFields() + logEvent("Initial Adapt 2 before Solve") + self.PUMI_adaptMesh() + #NS_base has a fairly complicated time stepping loop structure #to accommodate fairly general split operator approaches. The #outer loop is over user defined time intervals for the entire @@ -1430,7 +1569,7 @@ def calculateSolution(self,runName): self.nSequenceSteps = 0 nSequenceStepsLast=self.nSequenceSteps # prevent archiving the same solution twice - self.nSolveSteps=self.nList[0].adaptMesh_nSteps-1 + self.nSolveSteps=0 for (self.tn_last,self.tn) in zip(self.tnList[:-1],self.tnList[1:]): logEvent("==============================================================",level=0) logEvent("Solving over interval [%12.5e,%12.5e]" % (self.tn_last,self.tn),level=0) @@ -1460,9 +1599,11 @@ def calculateSolution(self,runName): lm.u[ci].dof_last[:] = lm.u[ci].dof logEvent("saving previous velocity dofs %s" % self.nSolveSteps) + logEvent("Split operator iteration %i" % (self.systemStepController.its,),level=3) self.nSequenceSteps += 1 for (self.t_stepSequence,model) in self.systemStepController.stepSequence: + logEvent("NumericalAnalytics Model %s " % (model.name), level=0) logEvent("Model: %s" % (model.name),level=1) logEvent("NumericalAnalytics Time Step " + `self.t_stepSequence`, level=0) @@ -1472,6 +1613,7 @@ def calculateSolution(self,runName): m.t_mesh = self.systemStepController.t_system_last m.updateAfterMeshMotion() m.tLast_mesh = m.t_mesh + self.preStep(model) self.setWeakDirichletConditions(model) @@ -1490,13 +1632,14 @@ def calculateSolution(self,runName): logEvent("Model substep t=%12.5e for model %s" % (self.tSubstep,model.name),level=3) #TODO: model.stepController.substeps doesn't seem to be updated after a solver failure unless model.stepController.stepExact is true logEvent("Model substep t=%12.5e for model %s model.timeIntegration.t= %12.5e" % (self.tSubstep,model.name,model.levelModelList[-1].timeIntegration.t),level=3) + model.stepController.setInitialGuess(model.uList,model.rList) - solverFailed = model.solver.solveMultilevel(uList=model.uList, rList=model.rList, par_uList=model.par_uList, par_rList=model.par_rList) + Profiling.memory("solver.solveMultilevel") if self.opts.wait: raw_input("Hit any key to continue") diff --git a/proteus/TimeIntegration.py b/proteus/TimeIntegration.py index 8b5fb45f35..d459b5a073 100644 --- a/proteus/TimeIntegration.py +++ b/proteus/TimeIntegration.py @@ -221,6 +221,12 @@ def setFromOptions(self,nOptions): allow classes to set various numerical parameters """ pass + def postAdaptUpdate(self,oldTime): + """ + ensure array histories are properly set after adapt. This should only work for backwards euler and vbdf at the moment. + """ + print "WARNING: calling abstract function that doesn't have a definition for your object so just passing for now" + pass NoIntegration = TI_base @@ -355,6 +361,14 @@ def setFromOptions(self,nOptions): """ if 'runCFL' in dir(nOptions): self.runCFL = nOptions.runCFL + def postAdaptUpdate(self,oldTime): + """This looks exactly like updateTimeHistory with the exception of setting self.tLast=self.t""" + for ci in self.massComponents: + self.m_last[ci].flat[:] = self.m_tmp[ci].flat + + #unclear why this is needed; a change of machine precision in these values yield larger deviations in the gauges + self.dtLast = oldTime.dtLast + class SSP(BackwardEuler_cfl): def __init__(self,transport,runCFL=0.9,integrateInterpolationPoints=False): @@ -1478,6 +1492,28 @@ def updateTimeHistory(self,resetFromDOF=False): #decide where this goes self.chooseOrder() # + def postAdaptUpdate(self,oldTime): + """This looks exactly like updateTimeHistory with the exeption of setting self.tLast=self.t""" + self.nUpdatesTimeHistoryCalled = oldTime.nUpdatesTimeHistoryCalled + for n in range(self.max_order-1): + for ci in self.massComponents: + self.m_history[n+1][ci].flat[:] =self.m_history[n][ci].flat + self.mt_history[n+1][ci].flat[:]=self.mt_history[n][ci].flat + # + self.dt_history[n+1] = self.dt_history[n] + # + self.dt_history[0] = self.dt + for ci in self.massComponents: + self.m_history[0][ci].flat[:] = self.m_tmp[ci].flat + self.mt_history[0][ci].flat[:]= self.mt_tmp[ci].flat + + self.needToCalculateBDFCoefs = True + #decide where this goes + self.chooseOrder() + + #unclear why this is needed; a change of machine precision in these values yield larger deviations in the gauges + self.dt_history[:] = oldTime.dt_history[:] + def computeErrorEstimate(self): """calculate :math:`\vec{e}^{n+1}` diff --git a/proteus/TransportCoefficients.py b/proteus/TransportCoefficients.py index 66063b2c3c..073c2e54c2 100644 --- a/proteus/TransportCoefficients.py +++ b/proteus/TransportCoefficients.py @@ -242,6 +242,11 @@ def postStep(self,t,firstStep=False): Give the TC object an opportunity to modify itself after the time step. """ return {} + def postAdaptStep(self): + """ + Give opportunity to TC object to maintain consistency with previous TC object + """ + return {} def allocateDummyCoefficients(self,c=None,nPoints=101,uMax=1.0,uMin=0.0,nSpace=1): """ Allocate some coefficient dictionaries to use for viewing the coefficients diff --git a/proteus/mprans/MCorr.py b/proteus/mprans/MCorr.py index 646096355e..1a48eb9e03 100644 --- a/proteus/mprans/MCorr.py +++ b/proteus/mprans/MCorr.py @@ -166,11 +166,14 @@ def preStep(self, t, firstStep=False): def postStep(self, t, firstStep=False): if self.applyCorrection: # ls + self.lsModel.u[0].dof += self.massCorrModel.u[0].dof self.lsModel.q[('u', 0)] += self.massCorrModel.q[('u', 0)] self.lsModel.ebqe[('u', 0)] += self.massCorrModel.ebqe[('u', 0)] self.lsModel.q[('grad(u)', 0)] += self.massCorrModel.q[('grad(u)', 0)] self.lsModel.ebqe[('grad(u)', 0)] += self.massCorrModel.ebqe[('grad(u)', 0)] + + # vof if self.edgeBasedStabilizationMethods == False: self.massCorrModel.setMassQuadrature() @@ -189,10 +192,16 @@ def postStep(self, t, firstStep=False): self.massCorrModel.mesh.nElements_owned),), level=2) logEvent("Phase 0 mass (consistent) after mass correction (LS) %12.5e" % (self.massCorrModel.calculateMass(self.lsModel.q[('m', 0)]),), level=2) copyInstructions = {} + # get the waterline on the obstacle if option set in NCLS (boundary==7) self.lsModel.computeWaterline(t) return copyInstructions + def postAdaptStep(self): + if self.applyCorrection: + # ls + self.lsModel.ebqe[('grad(u)', 0)][:] = self.massCorrModel.ebqe[('grad(u)', 0)] + def evaluate(self, t, c): import math if c[('u', 0)].shape == self.q_u_ls.shape: @@ -672,6 +681,7 @@ def getResidual(self, u, r): self.mesh.exteriorElementBoundariesArray, self.mesh.elementBoundaryElementsArray, self.mesh.elementBoundaryLocalElementBoundariesArray) + logEvent("Global residual", level=9, data=r) self.coefficients.massConservationError = fabs(globalSum(r[:self.mesh.nNodes_owned].sum())) logEvent(" Mass Conservation Error: ", level=3, data=self.coefficients.massConservationError) @@ -1137,6 +1147,7 @@ def setMassQuadratureEdgeBasedStabilizationMethods(self): self.lumped_L2p_vof_mass_correction.size) def setMassQuadrature(self): + self.mcorr.setMassQuadrature( # element self.u[0].femSpace.elementMaps.psi, self.u[0].femSpace.elementMaps.grad_psi, diff --git a/proteus/mprans/RDLS.py b/proteus/mprans/RDLS.py index f876fd7ac0..a7381cbb7c 100644 --- a/proteus/mprans/RDLS.py +++ b/proteus/mprans/RDLS.py @@ -1002,7 +1002,6 @@ def getResidual(self, u, r): if self.interface_locator[gi] == 1.0: self.u[0].dof[gi] = self.coefficients.dof_u0[gi] # END OF FREEZING INTERFACE # - self.calculateResidual( # element self.u[0].femSpace.elementMaps.psi, self.u[0].femSpace.elementMaps.grad_psi, @@ -1073,6 +1072,7 @@ def getResidual(self, u, r): self.lumped_qz, self.coefficients.alpha/self.elementDiameter.min()) + # FREEZE INTERFACE # if self.coefficients.freeze_interface_within_elliptic_redist==True: for gi in range(len(self.u[0].dof)): @@ -1241,12 +1241,6 @@ def estimate_mt(self): def calculateSolutionAtQuadrature(self): pass - # self.q[('u',0)].flat[:]=0.0 - # for eN in range(self.u[0].femSpace.elementMaps.mesh.nElements_global): - # for k in range(self.nQuadraturePoints_element): - # for j in self.u[0].femSpace.referenceFiniteElement.localFunctionSpace.range_dim: - # J = self.u[0].femSpace.dofMap.l2g[eN,j] - # self.q[('u',0)][eN,k]+=self.u[0].dof[J]*self.u[0].femSpace.psi[k,j] def calculateAuxiliaryQuantitiesAfterStep(self): pass diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/README b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/README new file mode 100644 index 0000000000..28dadeaa57 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/README @@ -0,0 +1,11 @@ +This test is to ensure that the restart process after adapt is working properly with + +1) RANS2P lagged subgrid error +2) RANS2P lagged shock capturing +3) lagged eddy-viscosity +4) Lagged shock capturing for VF and LS +5) LAG_LES=1.0 for dynamic smagorinsky + +This is for the nonconservative form of NS and for adaptive time stepping. +Oddly, fixed time stepping yields slight differences on the order of 1e-14 which throws off the current (simple) method for comparison. + diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/Reconstructed.dmg b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/Reconstructed.dmg new file mode 100644 index 0000000000..427ff11a38 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/Reconstructed.dmg @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2086854413c5dc58922f94e37197d2a52939765489d26013ab2bb63074516c95 +size 98 diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/dambreak_Colagrossi.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/dambreak_Colagrossi.py new file mode 100644 index 0000000000..22f9c27118 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/dambreak_Colagrossi.py @@ -0,0 +1,311 @@ +""" +Dambreak flow - Collagrosi and Landrini (2003) +""" +import numpy as np +from math import sqrt +from proteus.default_n import * +from proteus import (Domain, Context, + FemTools as ft, + #SpatialTools as st, + MeshTools as mt, + WaveTools as wt, + Gauges) +from proteus.mprans import SpatialTools as st +from proteus.Profiling import logEvent +from proteus.mprans.SpatialTools import Tank2D +from proteus.Gauges import PointGauges + + +# predefined options +opts=Context.Options([ + # water column + ("water_level", 0.6, "Height of water column in m"), + ("water_width", 1.2, "Width of water column in m"), + # tank + ("tank_dim", (1.8, 1.0), "Dimensions of the tank in m"), + #gravity + ("g",(0,-9.81,0), "Gravity vector in m/s^2"), + # gauges + ("gauge_output", True, "Produce gauge data"), + ("gauge_location_p", (1.8, 0.1, 0), "Pressure gauge location in m"), + # mesh refinement and timestep + ("refinement", 32 ,"Refinement level, he = L/(4*refinement - 1), where L is the horizontal dimension"), + ("cfl", 0.33 ,"Target cfl"), + # run time options + #("T", 0.011,"Simulation time in s"), + ("T", 0.05,"Simulation time in s"), + ("dt_fixed", 0.01, "Fixed time step in s"), + ("dt_init", 0.001 ,"Maximum initial time step in s"), + ("useHex", False, "Use a hexahedral structured mesh"), + ("structured", False, "Use a structured triangular mesh"), + ("gen_mesh", False ,"Generate new mesh"), + ("usePUMI", False ,"Generate new mesh"), + ("adapt",0,"adapt"), + ("fixedTimeStep",True,"choice of time stepping") + ]) + + +pressure_gauges = PointGauges(gauges=((('p',), + ((1.8,0.1,0.0), + (1.4,0.0,0.0) + )),), + fileName="pressure.csv") + + +# ----- CONTEXT ------ # + +# water +waterLine_z = opts.water_level +waterLine_x = opts.water_width + +# tank +tank_dim = opts.tank_dim + +########################################## +# Discretization Input Options # +########################################## + +#[temp] temporary location +backgroundDiffusionFactor = 0.01 + +refinement = opts.refinement +genMesh = opts.gen_mesh +usePUMI = opts.usePUMI +movingDomain = False +checkMass = False +useOldPETSc = False +useSuperlu = False +timeDiscretization = 'be' # 'vbdf', 'be', 'flcbdf' +#timeDiscretization = 'vbdf' # 'vbdf', 'be', 'flcbdf' +spaceOrder = 1 +useHex = opts.useHex +structured = opts.structured +useRBLES = 0.0 +useMetrics = 1.0 +applyRedistancing = False#True +applyCorrection = False#True +useVF = 0.0 +useOnlyVF = False +useRANS = 0 # 0 -- None + # 1 -- K-Epsilon + # 2 -- K-Omega + +# ----- INPUT CHECKS ----- # +if spaceOrder not in [1,2]: + raise ValueError("INVALID: spaceOrder(" + str(spaceOrder) + ")") + +if useRBLES not in [0.0, 1.0]: + raise ValueError("INVALID: useRBLES(" + str(useRBLES) + ")") + +if useMetrics not in [0.0, 1.0]: + raise ValueError("INVALID: useMetrics(" + str(useMetrics) + ")") + +# ----- DISCRETIZATION ----- # +nd = 2 +if spaceOrder == 1: + hFactor = 1.0 + if useHex: + basis = ft.C0_AffineLinearOnCubeWithNodalBasis + elementQuadrature = ft.CubeGaussQuadrature(nd, 2) + elementBoundaryQuadrature = ft.CubeGaussQuadrature(nd - 1, 2) + else: + basis = ft.C0_AffineLinearOnSimplexWithNodalBasis + elementQuadrature = ft.SimplexGaussQuadrature(nd, 3) + elementBoundaryQuadrature = ft.SimplexGaussQuadrature(nd - 1, 3) +elif spaceOrder == 2: + hFactor = 0.5 + if useHex: + basis = ft.C0_AffineLagrangeOnCubeWithNodalBasis + elementQuadrature = ft.CubeGaussQuadrature(nd, 4) + elementBoundaryQuadrature = ft.CubeGaussQuadrature(nd - 1, 4) + else: + basis = ft.C0_AffineQuadraticOnSimplexWithNodalBasis + elementQuadrature = ft.SimplexGaussQuadrature(nd, 4) + elementBoundaryQuadrature = ft.SimplexGaussQuadrature(nd - 1, 4) + +########################################## +# Numerical Options and Other Parameters # +########################################## + +weak_bc_penalty_constant = 100.0 +nLevels = 1 + +# ----- PHYSICAL PROPERTIES ----- # + +# Water +rho_0 = 998.2 +nu_0 = 1.004e-6 + +# Air +rho_1 = 1.205 +nu_1 = 1.500e-5 + +# Surface Tension +sigma_01 = 0.0 + +# Gravity +g = opts.g + +# ----- TIME STEPPING & VELOCITY----- # + +T = opts.T +dt_fixed = opts.dt_fixed +dt_init = min(0.1 * dt_fixed, opts.dt_init) +runCFL = opts.cfl +nDTout = int(round(T / dt_fixed)) +fixedTimeStep = opts.fixedTimeStep + +# ----- DOMAIN ----- # + +# domain = Domain.PlanarStraightLineGraphDomain() +# domain replacement +if useHex: + nnx = 4 * refinement + 1 + nny = 2 * refinement + 1 + hex = True + domain = Domain.RectangularDomain(tank_dim) +elif structured: + nnx = 4 * refinement + nny = 2 * refinement + domain = Domain.RectangularDomain(tank_dim) + boundaryTags = domain.boundaryTags + +elif usePUMI and not genMesh: + from proteus.MeshAdaptPUMI import MeshAdaptPUMI + domain = Domain.PUMIDomain(dim=nd) #initialize the domain + #boundaryTags=baseDomain.boundaryFlags + he = 0.06 + adaptMeshFlag = opts.adapt#1 + adaptMesh_nSteps =3#5 + adaptMesh_numIter = 5 + hmax = he; + hmin = he/4.0; + hPhi = he/2.0;#/4.0 + domain.PUMIMesh=MeshAdaptPUMI.MeshAdaptPUMI(hmax=hmax, hmin=hmin, hPhi = hPhi, adaptMesh=adaptMeshFlag, numIter=adaptMesh_numIter, numAdaptSteps=adaptMesh_nSteps, sfConfig="pseudo",logType="off",reconstructedFlag=2,gradingFact=1.2) + domain.PUMIMesh.loadModelAndMesh("Reconstructed.dmg", "Reconstructed.smb") + +else: + domain = Domain.PlanarStraightLineGraphDomain() + +if genMesh and usePUMI: + from proteus.MeshAdaptPUMI import MeshAdaptPUMI + domain.PUMIMesh=MeshAdaptPUMI.MeshAdaptPUMI() + + +# ----- TANK ----- # + +tank = Tank2D(domain, tank_dim) + +# ----- EXTRA BOUNDARY CONDITIONS ----- # + +tank.BC['y+'].setAtmosphere() +tank.BC['y-'].setFreeSlip() +tank.BC['x+'].setFreeSlip() +tank.BC['x-'].setFreeSlip() + +# ----- MESH CONSTRUCTION ----- # + +he = 0.06#tank_dim[0] / float(4 * refinement - 1) +domain.MeshOptions.he = he +st.assembleDomain(domain) + +# ----- STRONG DIRICHLET ----- # + +ns_forceStrongDirichlet = False + +# ----- NUMERICAL PARAMETERS ----- # + +if useMetrics: + ns_shockCapturingFactor = 0.25 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.25 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.25 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.25 + rd_lag_shockCapturing = True + epsFact_density = epsFact_viscosity = epsFact_curvature \ + = epsFact_vof = ecH \ + = epsFact_consrv_dirac = epsFact_density \ + = 3.0 + epsFact_redistance = 0.33 + epsFact_consrv_diffusion = 0.1 + redist_Newton = True + kappa_shockCapturingFactor = 0.25 + kappa_lag_shockCapturing = False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.25 + dissipation_lag_shockCapturing = False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 +else: + ns_shockCapturingFactor = 0.9 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.9 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.9 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.9 + rd_lag_shockCapturing = False + epsFact_density = epsFact_viscosity = epsFact_curvature \ + = epsFact_vof = ecH \ + = epsFact_consrv_dirac = epsFact_density \ + = 1.5 + epsFact_redistance = 0.33 + epsFact_consrv_diffusion = 1.0 + redist_Newton = False + kappa_shockCapturingFactor = 0.9 + kappa_lag_shockCapturing = True #False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.9 + dissipation_lag_shockCapturing = True #False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 + +# ----- NUMERICS: TOLERANCES ----- # + +ns_nl_atol_res = max(1.0e-10, 0.001 * he ** 2) +vof_nl_atol_res = max(1.0e-10, 0.001 * he ** 2) +ls_nl_atol_res = max(1.0e-10, 0.001 * he ** 2) +rd_nl_atol_res = max(1.0e-10, 0.005 * he) +mcorr_nl_atol_res = max(1.0e-10, 0.001 * he ** 2) +kappa_nl_atol_res = max(1.0e-10, 0.001 * he ** 2) +dissipation_nl_atol_res = max(1.0e-10, 0.001 * he ** 2) + +# ----- TURBULENCE MODELS ----- # + +ns_closure = 2 #1-classic smagorinsky, 2-dynamic smagorinsky, 3 -- k-epsilon, 4 -- k-omega +if useRANS == 1: + ns_closure = 3 +elif useRANS == 2: + ns_closure = 4 + +########################################## +# Initial conditions for free-surface # +########################################## + +def signedDistance(x): + phi_x = x[0] - waterLine_x + phi_z = x[1] - waterLine_z + if phi_x < 0.0: + if phi_z < 0.0: + return max(phi_x, phi_z) + else: + return phi_z + else: + if phi_z < 0.0: + return phi_x + else: + return sqrt(phi_x ** 2 + phi_z ** 2) diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/dambreak_Colagrossi_so.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/dambreak_Colagrossi_so.py new file mode 100644 index 0000000000..ae729f1507 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/dambreak_Colagrossi_so.py @@ -0,0 +1,63 @@ +""" +Split operator module for two-phase flow +""" + +import os +from proteus.default_so import * +from proteus import Context + +# Create context from main module +name_so = os.path.basename(__file__) +if '_so.py' in name_so[-6:]: + name = name_so[:-6] +elif '_so.pyc' in name_so[-7:]: + name = name_so[:-7] +else: + raise NameError('Split operator module must end with "_so.py"') + +try: + case = __import__(name) + Context.setFromModule(case) + ct = Context.get() +except ImportError: + raise ImportError(str(name) + '.py not found') + + +# List of p/n files +pnList = [] + +if ct.useOnlyVF: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n")] +else: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n"), + ("ls_p", "ls_n"), + ("redist_p", "redist_n"), + ("ls_consrv_p", "ls_consrv_n") + ] + +if ct.useRANS > 0: + pnList.append(("kappa_p", + "kappa_n")) + pnList.append(("dissipation_p", + "dissipation_n")) +name = "dambreak_Colagrossi_p" + +if ct.timeDiscretization == 'flcbdf': + systemStepControllerType = Sequential_MinFLCBDFModelStep + systemStepControllerType = Sequential_MinAdaptiveModelStep +else: + systemStepControllerType = Sequential_MinAdaptiveModelStep + +#systemStepControllerType = Sequential_FixedStep #Sequential_FixedStep_Simple # uses time steps in so.tnList +if(ct.fixedTimeStep): + systemStepControllerType = Sequential_FixedStep_Simple # uses time steps in so.tnList +#dt_system_fixed = 0.01; +systemStepExact=False; + + +needEBQ_GLOBAL = False +needEBQ = False + +tnList=[0.0,ct.dt_init]+[ct.dt_init+ i*ct.dt_fixed for i in range(1,ct.nDTout+1)] diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_consrv_n.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_consrv_n.py new file mode 100644 index 0000000000..a57747485c --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_consrv_n.py @@ -0,0 +1,73 @@ +from proteus.default_n import * +from proteus import (StepControl, + TimeIntegration, + NonlinearSolvers, + LinearSolvers, + LinearAlgebraTools, + NumericalFlux) +import ls_consrv_p as physics +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + +#time stepping +runCFL = ct.runCFL +timeIntegrator = TimeIntegration.ForwardIntegrator +timeIntegration = TimeIntegration.NoIntegration + +#mesh options +nLevels = ct.nLevels +parallelPartitioningType = mesh.parallelPartitioningType +nLayersOfOverlapForParallel = mesh.nLayersOfOverlapForParallel +restrictFineSolutionToAllMeshes = mesh.restrictFineSolutionToAllMeshes +triangleOptions = mesh.triangleOptions + +elementQuadrature = ct.elementQuadrature +elementBoundaryQuadrature = ct.elementBoundaryQuadrature + +femSpaces = {0: ct.basis} + +subgridError = None +massLumping = False +numericalFluxType = NumericalFlux.DoNothing +conservativeFlux = None +shockCapturing = None + +fullNewtonFlag = True +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = LinearAlgebraTools.SparseMatrix + +if ct.useOldPETSc: + multilevelLinearSolver = LinearSolvers.PETSc + levelLinearSolver = LinearSolvers.PETSc +else: + multilevelLinearSolver = LinearSolvers.KSP_petsc4py + levelLinearSolver = LinearSolvers.KSP_petsc4py + +if ct.useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU + +linear_solver_options_prefix = 'mcorr_' +nonlinearSolverConvergenceTest = 'rits' +levelNonlinearSolverConvergenceTest = 'rits' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*ct.mcorr_nl_atol_res +nl_atol_res = ct.mcorr_nl_atol_res +useEisenstatWalker = False#True + +maxNonlinearIts = 50 +maxLineSearches = 0 + +auxiliaryVariables = ct.domain.auxiliaryVariables['ls_consrv'] diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_consrv_p.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_consrv_p.py new file mode 100644 index 0000000000..fc938a2e31 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_consrv_p.py @@ -0,0 +1,37 @@ +from proteus.default_p import * +from proteus.mprans import MCorr +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + + +genMesh = mesh.genMesh +movingDomain = ct.movingDomain +T = ct.T + +LevelModelType = MCorr.LevelModel + +coefficients = MCorr.Coefficients(LSModel_index=int(ct.movingDomain)+2, + V_model=int(ct.movingDomain)+0, + me_model=int(ct.movingDomain)+4, + VOFModel_index=int(ct.movingDomain)+1, + applyCorrection=ct.applyCorrection, + nd=nd, + checkMass=True, + useMetrics=ct.useMetrics, + epsFactHeaviside=ct.ecH, + epsFactDirac=ct.epsFact_consrv_dirac, + epsFactDiffusion=ct.epsFact_consrv_diffusion) + +class zero_phi: + def __init__(self): + pass + def uOfX(self, X): + return 0.0 + def uOfXT(self, X, t): + return 0.0 + +initialConditions = {0: zero_phi()} diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_n.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_n.py new file mode 100644 index 0000000000..57f554b2f0 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_n.py @@ -0,0 +1,92 @@ +from proteus.default_n import * +import ls_p as physics +from proteus import (StepControl, + TimeIntegration, + NonlinearSolvers, + LinearSolvers, + LinearAlgebraTools) +from proteus.mprans import NCLS +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + +if ct.timeDiscretization=='vbdf': + timeIntegration = TimeIntegration.VBDF + timeOrder=2 + stepController = StepControl.Min_dt_cfl_controller +elif ct.timeDiscretization=='flcbdf': + timeIntegration = TimeIntegration.FLCBDF + #stepController = FLCBDF_controller + stepController = StepControl.Min_dt_cfl_controller + time_tol = 10.0*ct.ls_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = TimeIntegration.BackwardEuler_cfl + stepController = StepControl.Min_dt_cfl_controller + +# time stepping +runCFL = ct.runCFL + +# mesh options +nLevels = ct.nLevels +parallelPartitioningType = mesh.parallelPartitioningType +nLayersOfOverlapForParallel = mesh.nLayersOfOverlapForParallel +restrictFineSolutionToAllMeshes = mesh.restrictFineSolutionToAllMeshes +triangleOptions = mesh.triangleOptions + + + +elementQuadrature = ct.elementQuadrature +elementBoundaryQuadrature = ct.elementBoundaryQuadrature + +femSpaces = {0: ct.basis} + +massLumping = False +conservativeFlux = None +numericalFluxType = NCLS.NumericalFlux + +subgridError = NCLS.SubgridError(coefficients=physics.coefficients, + nd=ct.domain.nd) +shockCapturing = NCLS.ShockCapturing(physics.coefficients, + ct.domain.nd, + shockCapturingFactor=ct.ls_shockCapturingFactor, + lag=ct.ls_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = LinearAlgebraTools.SparseMatrix + +if ct.useOldPETSc: + multilevelLinearSolver = LinearSolvers.PETSc + levelLinearSolver = LinearSolvers.PETSc +else: + multilevelLinearSolver = LinearSolvers.KSP_petsc4py + levelLinearSolver = LinearSolvers.KSP_petsc4py +if ct.useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU + +linear_solver_options_prefix = 'ncls_' +nonlinearSolverConvergenceTest = 'rits' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.0 +l_atol_res = 0.1*ct.ls_nl_atol_res +nl_atol_res = ct.ls_nl_atol_res +useEisenstatWalker = False#True + +maxNonlinearIts = 50 +maxLineSearches = 0 + +auxiliaryVariables = ct.domain.auxiliaryVariables['ls'] diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_p.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_p.py new file mode 100644 index 0000000000..8d4d520fe9 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/ls_p.py @@ -0,0 +1,35 @@ +from proteus.default_p import * +from proteus.mprans import NCLS +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = domain.nd +mesh = domain.MeshOptions + + +genMesh = mesh.genMesh +movingDomain = ct.movingDomain +T = ct.T + +LevelModelType = NCLS.LevelModel + +coefficients = NCLS.Coefficients(V_model=int(ct.movingDomain)+0, + RD_model=int(ct.movingDomain)+3, + ME_model=int(ct.movingDomain)+2, + checkMass=False, + useMetrics=ct.useMetrics, + epsFact=ct.ecH, + sc_uref=ct.ls_sc_uref, + sc_beta=ct.ls_sc_beta, + movingDomain=ct.movingDomain) + +dirichletConditions = {0: lambda x, flag: None} +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0: {}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return ct.signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/redist_n.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/redist_n.py new file mode 100644 index 0000000000..7a03ad77cb --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/redist_n.py @@ -0,0 +1,94 @@ +from proteus.default_n import * +from proteus import (StepControl, + TimeIntegration, + NonlinearSolvers, + LinearSolvers, + LinearAlgebraTools, + NumericalFlux) +from proteus.mprans import RDLS +import redist_p as physics +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + +# time stepping +runCFL = ct.runCFL + +# mesh options +nLevels = ct.nLevels +parallelPartitioningType = mesh.parallelPartitioningType +nLayersOfOverlapForParallel = mesh.nLayersOfOverlapForParallel +restrictFineSolutionToAllMeshes = mesh.restrictFineSolutionToAllMeshes +triangleOptions = mesh.triangleOptions + +elementQuadrature = ct.elementQuadrature +elementBoundaryQuadrature = ct.elementBoundaryQuadrature + +femSpaces = {0: ct.basis} + +massLumping = False +numericalFluxType = NumericalFlux.DoNothing +conservativeFlux = None +subgridError = RDLS.SubgridError(coefficients=physics.coefficients, + nd=ct.domain.nd) +shockCapturing = RDLS.ShockCapturing(coefficients=physics.coefficients, + nd=ct.domain.nd, + shockCapturingFactor=ct.rd_shockCapturingFactor, + lag=ct.rd_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +nonlinearSmoother = NonlinearSolvers.NLGaussSeidel +linearSmoother = None + +matrix = LinearAlgebraTools.SparseMatrix + +if ct.useOldPETSc: + multilevelLinearSolver = LinearSolvers.PETSc + levelLinearSolver = LinearSolvers.PETSc +else: + multilevelLinearSolver = LinearSolvers.KSP_petsc4py + levelLinearSolver = LinearSolvers.KSP_petsc4py +if ct.useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU + +if ct.redist_Newton: + timeIntegration = TimeIntegration.NoIntegration + stepController = StepControl.Newton_controller + maxNonlinearIts = 50 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'rits' + levelNonlinearSolverConvergenceTest = 'rits' + linearSolverConvergenceTest = 'r-true' + useEisenstatWalker = False +else: + timeIntegration = TimeIntegration.BackwardEuler_cfl + stepController = RDLS.PsiTC + runCFL=2.0 + psitc['nStepsForce']=3 + psitc['nStepsMax']=50 + psitc['reduceRatio']=10.0 + psitc['startRatio']=1.0 + rtol_res[0] = 0.0 + atol_res[0] = ct.rd_nl_atol_res + useEisenstatWalker = False#True + maxNonlinearIts = 1 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'rits' + levelNonlinearSolverConvergenceTest = 'rits' + linearSolverConvergenceTest = 'r-true' + +linear_solver_options_prefix = 'rdls_' + +nl_atol_res = ct.rd_nl_atol_res +tolFac = 0.0 +linTolFac = 0.0 +l_atol_res = 0.01*ct.rd_nl_atol_res + +auxiliaryVariables = ct.domain.auxiliaryVariables['redist'] diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/redist_p.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/redist_p.py new file mode 100644 index 0000000000..e4fd3db978 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/redist_p.py @@ -0,0 +1,41 @@ +from proteus.default_p import * +from proteus.mprans import RDLS +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = domain.nd +mesh = domain.MeshOptions + + +genMesh = mesh.genMesh +movingDomain = ct.movingDomain +T = ct.T + +""" +The redistancing equation in the sloshbox test problem. +""" + +LevelModelType = RDLS.LevelModel + +coefficients = RDLS.Coefficients(applyRedistancing=ct.applyRedistancing, + epsFact=ct.epsFact_redistance, + nModelId=int(ct.movingDomain)+2, + rdModelId=int(ct.movingDomain)+3, + useMetrics=ct.useMetrics, + backgroundDiffusionFactor=ct.backgroundDiffusionFactor) + +def getDBC_rd(x, flag): + pass + +dirichletConditions = {0: getDBC_rd} +weakDirichletConditions = {0: RDLS.setZeroLSweakDirichletBCsSimple} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0: {}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return ct.signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/test_MeshAdaptRestart.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/test_MeshAdaptRestart.py new file mode 100644 index 0000000000..b216661984 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/test_MeshAdaptRestart.py @@ -0,0 +1,37 @@ +#!/usr/bin/env python + +from nose.tools import eq_ as eq +from nose.tools import ok_ as ok +import subprocess +import os +import pytest + +@pytest.mark.slowTest + +def test_MeshAdaptRestart_adaptiveTime_BackwardEuler_baseline(verbose=0): + """Get baseline data""" + currentPath = os.path.dirname(os.path.abspath(__file__)) + runCommand = "cd "+currentPath+"; parun -C \"gen_mesh=False usePUMI=True adapt=0 fixedTimeStep=False\" -D \"baseline\" dambreak_Colagrossi_so.py;" + subprocess.call(runCommand,shell=True) + assert(True) + + +def test_MeshAdaptRestart_adaptiveTime_BackwardEuler(verbose=0): + """Test restart workflow""" + currentPath = os.path.dirname(os.path.abspath(__file__)) + runCommand = "cd "+currentPath+"; parun -C \"gen_mesh=False usePUMI=True adapt=1 fixedTimeStep=False\" -D \"adapt_0\" dambreak_Colagrossi_so.py;" + subprocess.call(runCommand,shell=True) + assert(True) + +def test_MeshAdaptRestart(verbose=0): + #subprocess.call("diff normal adapt_0/pressure.csv baseline/pressure.csv") + currentPath = os.path.dirname(os.path.abspath(__file__)) + with open(currentPath+'/baseline/pressure.csv') as file1, open(currentPath+'/adapt_0/pressure.csv') as file2: + for line1, line2 in zip(file1, file2): + if(line1 != line2): + pytest.fail("pressure gauge values are not the same!\n") + +if __name__ == '__main__': + import nose + nose.main(defaultTest='test_MeshAdaptRestart:test_MeshAdaptRestart_adaptiveTime_BackwardEuler_baseline, test_MeshAdaptRestart:test_MeshAdaptRestart_adaptiveTime_BackwardEuler, test_MeshAdaptRestart:test_MeshAdaptRestart') + diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/twp_navier_stokes_n.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/twp_navier_stokes_n.py new file mode 100644 index 0000000000..dc05070fa5 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/twp_navier_stokes_n.py @@ -0,0 +1,107 @@ +from proteus import (StepControl, + TimeIntegration, + NonlinearSolvers, + LinearSolvers, + LinearAlgebraTools) +from proteus.default_n import * +import twp_navier_stokes_p as physics +from proteus.mprans import RANS2P +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + +if ct.useHex or ct.structured: + nnx = ct.nnx + nny = ct.nny + + if ct.useHex: + quad = True + + +#time stepping +runCFL = ct.runCFL +if ct.timeDiscretization=='vbdf': + timeIntegration = TimeIntegration.VBDF + timeOrder=2 + stepController = StepControl.Min_dt_cfl_controller +elif ct.timeDiscretization=='flcbdf': + timeIntegration = TimeIntegration.FLCBDF + #stepController = FLCBDF_controller_sys + stepController = StepControl.Min_dt_cfl_controller + time_tol = 10.0*ct.ns_nl_atol_res + atol_u = {1:time_tol,2:time_tol} + rtol_u = {1:time_tol,2:time_tol} +else: + timeIntegration = TimeIntegration.BackwardEuler_cfl + stepController = StepControl.Min_dt_cfl_controller + +# mesh options +nLevels = ct.nLevels +parallelPartitioningType = mesh.parallelPartitioningType +nLayersOfOverlapForParallel = mesh.nLayersOfOverlapForParallel +restrictFineSolutionToAllMeshes = mesh.restrictFineSolutionToAllMeshes +triangleOptions = mesh.triangleOptions + +elementQuadrature = ct.elementQuadrature +elementBoundaryQuadrature = ct.elementBoundaryQuadrature + +femSpaces = {0: ct.basis, + 1: ct.basis, + 2: ct.basis} + +massLumping = False + +numericalFluxType = RANS2P.NumericalFlux +#numericalFluxType = NumericalFlux.DoNothing + +subgridError = RANS2P.SubgridError(coefficients=physics.coefficients, + nd=nd, + lag=ct.ns_lag_subgridError, + hFactor=ct.hFactor) +shockCapturing = RANS2P.ShockCapturing(coefficients=physics.coefficients, + nd=nd, + shockCapturingFactor=ct.ns_shockCapturingFactor, + lag=ct.ns_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +nonlinearSmoother = None +linearSmoother = LinearSolvers.SimpleNavierStokes2D + +matrix = LinearAlgebraTools.SparseMatrix + +if ct.useOldPETSc: + multilevelLinearSolver = LinearSolvers.PETSc + levelLinearSolver = LinearSolvers.PETSc +else: + multilevelLinearSolver = LinearSolvers.KSP_petsc4py + levelLinearSolver = LinearSolvers.KSP_petsc4py + +if ct.useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU + +linear_solver_options_prefix = 'rans2p_' +nonlinearSolverConvergenceTest = 'rits' +levelNonlinearSolverConvergenceTest = 'rits' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*ct.ns_nl_atol_res +nl_atol_res = ct.ns_nl_atol_res +useEisenstatWalker = False#True +maxNonlinearIts = 50 +maxLineSearches = 0 +if ct.useHex: + conservativeFlux = None +else: + conservativeFlux = None#{0: 'pwl-bdm-opt'} + +#auxiliaryVariables = ct.domain.auxiliaryVariables['twp'] +auxiliaryVariables=[ct.pressure_gauges] diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/twp_navier_stokes_p.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/twp_navier_stokes_p.py new file mode 100644 index 0000000000..b0c23a1e0c --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/twp_navier_stokes_p.py @@ -0,0 +1,88 @@ +from proteus.default_p import * +from proteus.mprans import RANS2P +import numpy as np +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = domain.nd +mesh = domain.MeshOptions + +genMesh = mesh.genMesh +movingDomain = ct.movingDomain +T = ct.T # might not be necessaryd + +LevelModelType = RANS2P.LevelModel +if ct.useOnlyVF: + LS_model = None +else: + LS_model = 2 +if ct.useRANS >= 1: + Closure_0_model = 5 + Closure_1_model = 6 + if ct.useOnlyVF: + Closure_0_model = 2 + Closure_1_model = 3 + if ct.movingDomain: + Closure_0_model += 1 + Closure_1_model += 1 +else: + Closure_0_model = None + Closure_1_model = None + +coefficients = RANS2P.Coefficients(epsFact=ct.epsFact_viscosity, + sigma=0.0, + rho_0=ct.rho_0, + nu_0=ct.nu_0, + rho_1=ct.rho_1, + nu_1=ct.nu_1, + g=ct.g, + nd=nd, + ME_model=int(ct.movingDomain)+0, + VF_model=int(ct.movingDomain)+1, + LS_model=int(ct.movingDomain)+LS_model, + Closure_0_model=Closure_0_model, + Closure_1_model=Closure_1_model, + epsFact_density=ct.epsFact_density, + stokes=False, + useVF=ct.useVF, + useRBLES=ct.useRBLES, + useMetrics=ct.useMetrics, + eb_adjoint_sigma=1.0, + eb_penalty_constant=ct.weak_bc_penalty_constant, + forceStrongDirichlet=ct.ns_forceStrongDirichlet, + turbulenceClosureModel=ct.ns_closure, + movingDomain=ct.movingDomain, + LAG_LES=1.0, + NONCONSERVATIVE_FORM=0.0) + +dirichletConditions = {0: lambda x, flag: domain.bc[flag].p_dirichlet.init_cython(), + 1: lambda x, flag: domain.bc[flag].u_dirichlet.init_cython(), + 2: lambda x, flag: domain.bc[flag].v_dirichlet.init_cython()} + +advectiveFluxBoundaryConditions = {0: lambda x, flag: domain.bc[flag].p_advective.init_cython(), + 1: lambda x, flag: domain.bc[flag].u_advective.init_cython(), + 2: lambda x, flag: domain.bc[flag].v_advective.init_cython()} + +diffusiveFluxBoundaryConditions = {0: {}, + 1: {1: lambda x, flag: domain.bc[flag].u_diffusive.init_cython()}, + 2: {2: lambda x, flag: domain.bc[flag].v_diffusive.init_cython()}} + +class PerturbedSurface_p: + def __init__(self,waterLevel): + self.waterLevel=waterLevel + def uOfXT(self,x,t): + if ct.signedDistance(x) < 0: + return -(L[1] - self.waterLevel)*ct.rho_1*ct.g[1] - (self.waterLevel - x[1])*ct.rho_0*ct.g[1] + else: + return -(L[1] - self.waterLevel)*ct.rho_1*ct.g[1] + +class AtRest: + def __init__(self): + pass + def uOfXT(self,x,t): + return 0.0 + +initialConditions = {0:PerturbedSurface_p(ct.waterLine_z), + 1:AtRest(), + 2:AtRest()} diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/vof_n.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/vof_n.py new file mode 100644 index 0000000000..62fcbd5b64 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/vof_n.py @@ -0,0 +1,94 @@ +from proteus.default_n import * +from proteus import (StepControl, + TimeIntegration, + NonlinearSolvers, + LinearSolvers, + LinearAlgebraTools) +import vof_p as physics +from proteus.mprans import VOF +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + +if ct.useHex or ct.structured: + nnx = ct.nnx + nny = ct.nny + +# time stepping +runCFL = ct.runCFL +if ct.timeDiscretization=='vbdf': + timeIntegration = TimeIntegration.VBDF + timeOrder=2 + stepController = StepControl.Min_dt_cfl_controller +elif ct.timeDiscretization=='flcbdf': + timeIntegration = TimeIntegration.FLCBDF + #stepController = FLCBDF_controller + stepController = StepControl.Min_dt_cfl_controller + time_tol = 10.0*ct.vof_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = TimeIntegration.BackwardEuler_cfl + stepController = StepControl.Min_dt_cfl_controller + +# mesh options +nLevels = ct.nLevels +parallelPartitioningType = mesh.parallelPartitioningType +nLayersOfOverlapForParallel = mesh.nLayersOfOverlapForParallel +restrictFineSolutionToAllMeshes = mesh.restrictFineSolutionToAllMeshes +triangleOptions = mesh.triangleOptions + +elementQuadrature = ct.elementQuadrature +elementBoundaryQuadrature = ct.elementBoundaryQuadrature + +femSpaces = {0: ct.basis} + +massLumping = False +numericalFluxType = VOF.NumericalFlux +conservativeFlux = None +subgridError = VOF.SubgridError(coefficients=physics.coefficients, + nd=ct.domain.nd) +shockCapturing = VOF.ShockCapturing(coefficients=physics.coefficients, + nd=ct.domain.nd, + shockCapturingFactor=ct.vof_shockCapturingFactor, + lag=ct.vof_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = LinearAlgebraTools.SparseMatrix + +if ct.useOldPETSc: + multilevelLinearSolver = LinearSolvers.PETSc + levelLinearSolver = LinearSolvers.PETSc +else: + multilevelLinearSolver = LinearSolvers.KSP_petsc4py + levelLinearSolver = LinearSolvers.KSP_petsc4py + +if ct.useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU + +linear_solver_options_prefix = 'vof_' +nonlinearSolverConvergenceTest = 'rits' +levelNonlinearSolverConvergenceTest = 'rits' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.0 +l_atol_res = 0.1*ct.vof_nl_atol_res +nl_atol_res = ct.vof_nl_atol_res +useEisenstatWalker = False#True + +maxNonlinearIts = 50 +maxLineSearches = 0 + +auxiliaryVariables = ct.domain.auxiliaryVariables['vof'] + diff --git a/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/vof_p.py b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/vof_p.py new file mode 100644 index 0000000000..f0ceb76020 --- /dev/null +++ b/proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/vof_p.py @@ -0,0 +1,45 @@ +from proteus.default_p import * +from proteus.ctransportCoefficients import smoothedHeaviside +from proteus.mprans import VOF +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = domain.nd +mesh = domain.MeshOptions + + +genMesh = mesh.genMesh +movingDomain = ct.movingDomain +T = ct.T + +LevelModelType = VOF.LevelModel +if ct.useOnlyVF: + RD_model = None + LS_model = None +else: + RD_model = 3 + LS_model = 2 + +coefficients = VOF.Coefficients(#LS_model=int(ct.movingDomain)+LS_model, + V_model=int(ct.movingDomain)+0, + RD_model=int(ct.movingDomain)+RD_model, + ME_model=int(ct.movingDomain)+1, + checkMass=True, + useMetrics=ct.useMetrics, + epsFact=ct.epsFact_vof, + sc_uref=ct.vof_sc_uref, + sc_beta=ct.vof_sc_beta, + movingDomain=ct.movingDomain) + +dirichletConditions = {0: lambda x, flag: domain.bc[flag].vof_dirichlet.init_cython()} + +advectiveFluxBoundaryConditions = {0: lambda x, flag: domain.bc[flag].vof_advective.init_cython()} + +diffusiveFluxBoundaryConditions = {0: {}} + +class PerturbedSurface_H: + def uOfXT(self,x,t): + return smoothedHeaviside(ct.ecH*ct.he,ct.signedDistance(x)) + +initialConditions = {0:PerturbedSurface_H()} diff --git a/proteus/tests/MeshAdaptPUMI/test_gmshCheck.py b/proteus/tests/MeshAdaptPUMI/test_gmshCheck.py index b47a3f6702..f637991795 100755 --- a/proteus/tests/MeshAdaptPUMI/test_gmshCheck.py +++ b/proteus/tests/MeshAdaptPUMI/test_gmshCheck.py @@ -34,7 +34,8 @@ def test_gmshLoadAndAdapt(verbose=0): nu = numpy.array([1.004e-6, 1.004e-6]) g = numpy.asarray([0.0,0.0,0.0]) deltaT = 1.0 #dummy number - domain.PUMIMesh.transferPropertiesToPUMI(rho,nu,g,deltaT) + epsFact = 1.0 #dummy number + domain.PUMIMesh.transferPropertiesToPUMI(rho,nu,g,deltaT,epsFact) #Couette Flow Lz = 0.05 @@ -99,7 +100,8 @@ def test_2DgmshLoadAndAdapt(verbose=0): nu = numpy.array([1.004e-6, 1.004e-6]) g = numpy.asarray([0.0,0.0]) deltaT = 1.0 #dummy number - domain.PUMIMesh.transferPropertiesToPUMI(rho,nu,g,deltaT) + epsFact = 1.0 #dummy number + domain.PUMIMesh.transferPropertiesToPUMI(rho,nu,g,deltaT,epsFact) #Couette Flow Lz = 0.05 diff --git a/proteus/tests/MeshAdaptPUMI/test_poiseuilleError.py b/proteus/tests/MeshAdaptPUMI/test_poiseuilleError.py index d0ff000519..9f8c891ab1 100755 --- a/proteus/tests/MeshAdaptPUMI/test_poiseuilleError.py +++ b/proteus/tests/MeshAdaptPUMI/test_poiseuilleError.py @@ -38,7 +38,8 @@ def test_poiseuilleError(verbose=0): nu = numpy.array([1.004e-6, 1.004e-6]) g = numpy.asarray([0.0,0.0,0.0]) deltaT = 1.0 #dummy number - domain.PUMIMesh.transferPropertiesToPUMI(rho,nu,g,deltaT) + epsFact = 1.0 #dummy number + domain.PUMIMesh.transferPropertiesToPUMI(rho,nu,g,deltaT,epsFact) #Poiseuille Flow Ly=0.2 diff --git a/setup.py b/setup.py index 7134414ace..e0e6cda042 100644 --- a/setup.py +++ b/setup.py @@ -569,6 +569,7 @@ def setup_given_extensions(extensions): 'proteus.tests.matrix_constructor.import_modules', 'proteus.MeshAdaptPUMI', 'proteus.tests.MeshAdaptPUMI', + 'proteus.tests.MeshAdaptPUMI.gauge_compare.dambreak_Colagrossi_2D', 'proteus.tests.mesh_tests', 'proteus.tests.mesh_tests.import_modules', 'proteus.tests.poisson_2d', @@ -717,6 +718,9 @@ def setup_given_extensions(extensions): 'proteus/tests/MeshAdaptPUMI/Rectangle.dmg', 'proteus/tests/MeshAdaptPUMI/TwoQuads0.smb', 'proteus/tests/MeshAdaptPUMI/TwoQuads.dmg']), + (os.path.join(proteus_install_path,'tests','MeshAdaptPUMI','gauge_compare','dambreak_Colagrossi_2D'), + ['proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/Reconstructed.dmg', + 'proteus/tests/MeshAdaptPUMI/gauge_compare/dambreak_Colagrossi_2D/Reconstructed0.smb']), (os.path.join(proteus_install_path,'tests','poisson_2d'), ['proteus/tests/poisson_2d/square4x4.3dm', 'proteus/tests/poisson_2d/square4x4.bc']),