[cig-commits] commit: Compute maxDiffusivity once at the beginning from material properties. Also robustly handle the case where the diffusivity is 0
Mercurial
hg at geodynamics.org
Mon Nov 21 12:11:52 PST 2011
changeset: 817:ee8331550da9
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Mon Nov 21 12:11:43 2011 -0800
files: SLE/ProvidedSystems/AdvectionDiffusion/src/AdvectionDiffusionSLE.cxx SLE/ProvidedSystems/AdvectionDiffusion/src/AdvectionDiffusionSLE.h SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.cxx SLE/ProvidedSystems/AdvectionDiffusion/src/Timestep.cxx SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.cxx
description:
Compute maxDiffusivity once at the beginning from material properties. Also robustly handle the case where the diffusivity is 0
diff -r a4533511ac63 -r ee8331550da9 SLE/ProvidedSystems/AdvectionDiffusion/src/AdvectionDiffusionSLE.cxx
--- a/SLE/ProvidedSystems/AdvectionDiffusion/src/AdvectionDiffusionSLE.cxx Sun Nov 20 23:08:07 2011 -0800
+++ b/SLE/ProvidedSystems/AdvectionDiffusion/src/AdvectionDiffusionSLE.cxx Mon Nov 21 12:11:43 2011 -0800
@@ -427,8 +427,6 @@ void _AdvectionDiffusionSLE_Initialise(
/* Solution Vectors */
Stg_Component_Initialise( self->phiVector, data, False );
Stg_Component_Initialise( self->phiDotVector, data, False );
-
- AdvectionDiffusionSLE_ResetStoredValues( self );
}
void _AdvectionDiffusionSLE_Execute(void* sle, void* _context)
@@ -439,9 +437,7 @@ void _AdvectionDiffusionSLE_Execute(void
if(context->timeStep!=context->restartTimestep)
{
- AdvectionDiffusionSLE_ResetStoredValues(self);
self->currentDt=dt;
-
_SystemLinearEquations_Execute(self,context);
}
}
@@ -451,14 +447,3 @@ Vec _AdvectionDiffusionSLE_GetResidual(
AdvectionDiffusionSLE* self = (AdvectionDiffusionSLE*) sle;
return self->residual->vector;
}
-
-#define SMALL_VALUE 1.0e-99
-
-void AdvectionDiffusionSLE_ResetStoredValues( void* sle ) {
- AdvectionDiffusionSLE* self = (AdvectionDiffusionSLE*) sle;
-
- self->maxDiffusivity = SMALL_VALUE;
-}
-
-
-
diff -r a4533511ac63 -r ee8331550da9 SLE/ProvidedSystems/AdvectionDiffusion/src/AdvectionDiffusionSLE.h
--- a/SLE/ProvidedSystems/AdvectionDiffusion/src/AdvectionDiffusionSLE.h Sun Nov 20 23:08:07 2011 -0800
+++ b/SLE/ProvidedSystems/AdvectionDiffusion/src/AdvectionDiffusionSLE.h Mon Nov 21 12:11:43 2011 -0800
@@ -132,7 +132,5 @@
//Vector* _AdvectionDiffusionSLE_GetResidual( void* sle, Index fv_I );
Vec _AdvectionDiffusionSLE_GetResidual( void* sle, Index fv_I );
- void AdvectionDiffusionSLE_ResetStoredValues( void* sle );
-
#endif
diff -r a4533511ac63 -r ee8331550da9 SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.cxx
--- a/SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.cxx Sun Nov 20 23:08:07 2011 -0800
+++ b/SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.cxx Mon Nov 21 12:11:43 2011 -0800
@@ -56,6 +56,7 @@
#include <assert.h>
#include <string.h>
#include <stddef.h>
+#include <algorithm>
/* Textual name of this class */
const Type AdvDiffResidualForceTerm_Type = "AdvDiffResidualForceTerm";
@@ -245,70 +246,87 @@ void _AdvDiffResidualForceTerm_AssignFro
upwindFuncType );
}
-void _AdvDiffResidualForceTerm_Build( void* residual, void* data ) {
- AdvDiffResidualForceTerm* self = (AdvDiffResidualForceTerm*)residual;
- AdvDiffResidualForceTerm_MaterialExt* materialExt;
- Material_Index material_I;
- Material* material;
- Materials_Register* materials_Register = (Materials_Register*)(self->materials_Register);
- IntegrationPointsSwarm* swarm = (IntegrationPointsSwarm*)self->integrationSwarm;
- MaterialPointsSwarm** materialSwarms;
- Index materialSwarm_I;
- Stg_ComponentFactory* cf;
- char* name;
+void _AdvDiffResidualForceTerm_Build(void* residual, void* data)
+{
+ AdvDiffResidualForceTerm* self=(AdvDiffResidualForceTerm*)residual;
+ AdvDiffResidualForceTerm_MaterialExt* materialExt;
+ Material_Index material_I;
+ Material* material;
+ Materials_Register* materials_Register=
+ (Materials_Register*)(self->materials_Register);
+ IntegrationPointsSwarm* swarm=(IntegrationPointsSwarm*)self->integrationSwarm;
+ MaterialPointsSwarm** materialSwarms;
+ Index materialSwarm_I;
+ Stg_ComponentFactory* cf;
+ char* name;
- if(self->picSwarm)
- {
- swarm=(IntegrationPointsSwarm*)self->picSwarm;
- }
- else if(Stg_Class_IsInstance(swarm->mapper,OneToManyMapper_Type))
- {
- swarm=((OneToManyMapper*)(swarm->mapper))->swarm;
- }
- else if(Stg_Class_IsInstance(swarm->mapper,NearestNeighborMapper_Type))
- {
- swarm=((NearestNeighborMapper*)(swarm->mapper))->swarm;
- }
+ if(self->picSwarm)
+ {
+ swarm=(IntegrationPointsSwarm*)self->picSwarm;
+ }
+ else if(Stg_Class_IsInstance(swarm->mapper,OneToManyMapper_Type))
+ {
+ swarm=((OneToManyMapper*)(swarm->mapper))->swarm;
+ }
+ else if(Stg_Class_IsInstance(swarm->mapper,NearestNeighborMapper_Type))
+ {
+ swarm=((NearestNeighborMapper*)(swarm->mapper))->swarm;
+ }
- cf = self->context->CF;
+ cf = self->context->CF;
- _ForceTerm_Build( self, data );
+ _ForceTerm_Build( self, data );
- Stg_Component_Build( self->velocityField, data, False );
+ Stg_Component_Build( self->velocityField, data, False );
- /* Sort out material extension stuff */
- self->materialExtHandle = Materials_Register_AddMaterialExtension(
- self->materials_Register,
- self->type,
- sizeof(AdvDiffResidualForceTerm_MaterialExt) );
- for ( material_I = 0 ; material_I < Materials_Register_GetCount( materials_Register ) ; material_I++) {
- material = Materials_Register_GetByIndex( materials_Register, material_I );
- materialExt = (AdvDiffResidualForceTerm_MaterialExt*)ExtensionManager_GetFunc( material->extensionMgr, material, self->materialExtHandle );
+ AdvectionDiffusionSLE*
+ sle=Stg_CheckType(self->extraInfo,AdvectionDiffusionSLE);
+
+ /* Sort out material extension stuff */
+ self->materialExtHandle=
+ Materials_Register_AddMaterialExtension
+ (self->materials_Register, self->type,
+ sizeof(AdvDiffResidualForceTerm_MaterialExt));
+ for(material_I=0; material_I<Materials_Register_GetCount(materials_Register);
+ material_I++)
+ {
+ material=Materials_Register_GetByIndex(materials_Register,material_I);
+ materialExt=(AdvDiffResidualForceTerm_MaterialExt*)
+ ExtensionManager_GetFunc(material->extensionMgr,material,
+ self->materialExtHandle);
+
+ materialExt->diffusivity=
+ Stg_ComponentFactory_GetDouble(cf,material->name,"diffusivity",
+ self->defaultDiffusivity);
+ sle->maxDiffusivity=std::max(sle->maxDiffusivity,materialExt->diffusivity);
+ }
+
+ /* Create Swarm Variables of each material swarm this ip swarm is
+ mapped against */
+ materialSwarms=
+ IntegrationPointMapper_GetMaterialPointsSwarms(swarm->mapper,
+ &(self->materialSwarmCount));
+ self->diffusivitySwarmVariables=
+ (void**)Memory_Alloc_Array(MaterialSwarmVariable*,self->materialSwarmCount,
+ "DiffusivityVariables");
+
+ for(materialSwarm_I=0; materialSwarm_I<self->materialSwarmCount;
+ ++materialSwarm_I)
+ {
+ name=Stg_Object_AppendSuffix(materialSwarms[materialSwarm_I],
+ (Name)"Diffusivity");
+ self->diffusivitySwarmVariables[materialSwarm_I]=
+ MaterialSwarmVariable_New(name, (AbstractContext*)self->context,
+ materialSwarms[materialSwarm_I], 1,
+ (Materials_Register*)self->materials_Register,
+ self->materialExtHandle,
+ GetOffsetOfMember(*materialExt,diffusivity));
+ Memory_Free(name);
- materialExt->diffusivity = Stg_ComponentFactory_GetDouble( cf, material->name,
- (Dictionary_Entry_Key)"diffusivity",
- self->defaultDiffusivity );
- }
-
- /* Create Swarm Variables of each material swarm this ip swarm is mapped against */
- materialSwarms = IntegrationPointMapper_GetMaterialPointsSwarms(swarm->mapper, &(self->materialSwarmCount) );
- self->diffusivitySwarmVariables = (void**)Memory_Alloc_Array( MaterialSwarmVariable*, self->materialSwarmCount, "DiffusivityVariables" );
-
- for ( materialSwarm_I = 0; materialSwarm_I < self->materialSwarmCount; ++materialSwarm_I ) {
- name = Stg_Object_AppendSuffix( materialSwarms[materialSwarm_I], (Name)"Diffusivity" );
- self->diffusivitySwarmVariables[materialSwarm_I] = MaterialSwarmVariable_New(
- name,
- (AbstractContext*)self->context,
- materialSwarms[materialSwarm_I],
- 1,
- (Materials_Register*)self->materials_Register,
- self->materialExtHandle,
- GetOffsetOfMember( *materialExt, diffusivity ) );
- Memory_Free( name );
-
- /* Build new Swarm Variables */
- Stg_Component_Build( self->diffusivitySwarmVariables[materialSwarm_I], data, False );
- }
+ /* Build new Swarm Variables */
+ Stg_Component_Build(self->diffusivitySwarmVariables[materialSwarm_I],
+ data,False);
+ }
}
void _AdvDiffResidualForceTerm_Initialise( void* residual, void* data ) {
diff -r a4533511ac63 -r ee8331550da9 SLE/ProvidedSystems/AdvectionDiffusion/src/Timestep.cxx
--- a/SLE/ProvidedSystems/AdvectionDiffusion/src/Timestep.cxx Sun Nov 20 23:08:07 2011 -0800
+++ b/SLE/ProvidedSystems/AdvectionDiffusion/src/Timestep.cxx Mon Nov 21 12:11:43 2011 -0800
@@ -53,49 +53,60 @@
#include <string.h>
-double AdvectionDiffusionSLE_CalculateDt( void* advectionDiffusionSLE, FiniteElementContext* context ) {
- AdvectionDiffusionSLE* self = (AdvectionDiffusionSLE*) advectionDiffusionSLE;
- double advectionTimestep;
- double diffusionTimestep;
- double advectionTimestepGlobal;
- double diffusionTimestepGlobal;
- double timestep;
+double AdvectionDiffusionSLE_CalculateDt(void* advectionDiffusionSLE,
+ FiniteElementContext* context)
+{
+ AdvectionDiffusionSLE* self=(AdvectionDiffusionSLE*)advectionDiffusionSLE;
+ double advectionTimestep;
+ double diffusionTimestep;
+ double advectionTimestepGlobal;
+ double diffusionTimestepGlobal;
+ double timestep;
- Journal_DPrintf( self->debug, "In func: %s\n", __func__ );
+ Journal_DPrintf(self->debug, "In func: %s\n", __func__);
- /* It would be useful to introduce a limit to the change of step in here ... to prevent timesteps
- from becoming arbitrarily large in a single step */
+ /* It would be useful to introduce a limit to the change of step in
+ here ... to prevent timesteps from becoming arbitrarily large in
+ a single step */
- /* Calculate Courant Number */
- advectionTimestep = AdvectionDiffusionSLE_AdvectiveTimestep( self );
- diffusionTimestep = AdvectionDiffusionSLE_DiffusiveTimestep( self );
+ /* Calculate Courant Number */
+ advectionTimestep = AdvectionDiffusionSLE_AdvectiveTimestep(self);
+ diffusionTimestep = AdvectionDiffusionSLE_DiffusiveTimestep(self);
- MPI_Allreduce( &advectionTimestep, &advectionTimestepGlobal, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
- MPI_Allreduce( &diffusionTimestep, &diffusionTimestepGlobal, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
+ MPI_Allreduce(&advectionTimestep,&advectionTimestepGlobal,1,MPI_DOUBLE,
+ MPI_MIN,MPI_COMM_WORLD);
+ MPI_Allreduce(&diffusionTimestep,&diffusionTimestepGlobal,1,MPI_DOUBLE,
+ MPI_MIN,MPI_COMM_WORLD);
- Journal_DPrintf( self->debug, "%s Dominating. - Advective Timestep = %g - Diffusive Timestep = %g\n",
- advectionTimestepGlobal < diffusionTimestepGlobal ? "Advection" : "Diffusion",
- advectionTimestepGlobal, diffusionTimestepGlobal);
+ Journal_DPrintf(self->debug,
+ "%s Dominating. - Advective Timestep = %g - "
+ "Diffusive Timestep = %g\n",
+ advectionTimestepGlobal < diffusionTimestepGlobal
+ ? "Advection" : "Diffusion",
+ advectionTimestepGlobal, diffusionTimestepGlobal);
- /* Calculate Time Step */
- timestep = MIN( advectionTimestepGlobal, diffusionTimestepGlobal );
-
- return timestep;
+ /* Calculate Time Step */
+ if(diffusionTimestep<0)
+ timestep=advectionTimestepGlobal;
+ else
+ timestep=std::min(advectionTimestepGlobal,diffusionTimestepGlobal);
+ return timestep;
}
-double AdvectionDiffusionSLE_DiffusiveTimestep( void* advectionDiffusionSLE ) {
- AdvectionDiffusionSLE* self = (AdvectionDiffusionSLE*) advectionDiffusionSLE;
- double minSeparation;
- double minSeparationEachDim[3];
+double AdvectionDiffusionSLE_DiffusiveTimestep(void* advectionDiffusionSLE)
+{
+ AdvectionDiffusionSLE* self=(AdvectionDiffusionSLE*) advectionDiffusionSLE;
+ double minSeparation;
+ double minSeparationEachDim[3];
- Journal_DPrintf( self->debug, "In func: %s\n", __func__ );
+ Journal_DPrintf(self->debug,"In func: %s\n",__func__);
- FeVariable_GetMinimumSeparation( self->phiField, &minSeparation, minSeparationEachDim );
-
- /* This is quite a conservative estimate */
-
- return self->courantFactor * minSeparation * minSeparation / self->maxDiffusivity;
+ if(self->maxDiffusivity<=0)
+ return -1;
+ FeVariable_GetMinimumSeparation(self->phiField,&minSeparation,
+ minSeparationEachDim);
+ return self->courantFactor*minSeparation*minSeparation/self->maxDiffusivity;
}
diff -r a4533511ac63 -r ee8331550da9 SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.cxx
--- a/SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.cxx Sun Nov 20 23:08:07 2011 -0800
+++ b/SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.cxx Mon Nov 21 12:11:43 2011 -0800
@@ -117,11 +117,7 @@ double AdvDiffResidualForceTerm_UpwindDi
weight+=particle->weight;
}
averageDiffusivity/=weight;
-
- if(sle->maxDiffusivity<averageDiffusivity)
- sle->maxDiffusivity=averageDiffusivity;
-
-
+
/* Change Diffusivity if it is too small */
if(averageDiffusivity<SUPG_MIN_DIFFUSIVITY)
averageDiffusivity = SUPG_MIN_DIFFUSIVITY;
More information about the CIG-COMMITS
mailing list