[cig-commits] r13089 - long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3
echoi at geodynamics.org
echoi at geodynamics.org
Wed Oct 15 12:17:45 PDT 2008
Author: echoi
Date: 2008-10-15 12:17:45 -0700 (Wed, 15 Oct 2008)
New Revision: 13089
Modified:
long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/Force.c
long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/InitialConditions.c
long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/Register.c
Log:
minor changes
Modified: long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/Force.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/Force.c 2008-10-15 19:11:30 UTC (rev 13088)
+++ long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/Force.c 2008-10-15 19:17:45 UTC (rev 13089)
@@ -64,7 +64,6 @@
{
Snac_Context* context = (Snac_Context*)_context;
MeshLayout* meshLayout = (MeshLayout*)context->meshLayout;
- BlockGeometry* geometry = (BlockGeometry*)meshLayout->elementLayout->geometry;
double normal[3];
Node_ElementIndex nodeElement_I, nodeElementCount;
const double factor4 = 1.0f / 4.0f;
@@ -91,10 +90,13 @@
Material_Index material_I = element->material_I;
Snac_Material* material = &context->materialProperty[material_I];
Density phsDensity = material->phsDensity; /* node->density */
+ Density mantleDensity = 3300.0f;
double alpha = material->alpha;
double beta = material->beta;
- double drosubg = 500.0 * context->gravity;
+ double drosub = 0.0f;
+ double p_est = context->pisos + 0.5f * ( phsDensity + drosub ) * context->gravity * ( (*coord)[1] - element->rzbo );
+ double rosubg = context->gravity * ( phsDensity + drosub ) * ( 1.0 - alpha * (nodeT-material->reftemp) + beta * p_est );
double press_norm = 0.0f;
double area=0.0f, dhE=0.0f;
Normal* normal1;
@@ -141,7 +143,7 @@
/* Adjust all the signs to be consistent with this principle. */
/* Let's make isostatic pressure positive and normals to point to the direction of action, +y. */
/* dP should then have the same sign with dh. */
- press_norm = (float) ( drosubg * ( geometry->min[1] - dhE ) );
+ press_norm = context->pisos + rosubg * ( element->rzbo - dhE );
(*force)[0] += factor4 * ( press_norm * area * normal[0] );
(*force)[1] += factor4 * ( press_norm * area * normal[1] );
(*force)[2] += factor4 * ( press_norm * area * normal[2] );
@@ -503,10 +505,13 @@
Material_Index material_I = element->material_I;
Snac_Material* material = &context->materialProperty[material_I];
Density phsDensity = material->phsDensity; /* node->density */
+ Density mantleDensity = 3300.0f;
double alpha = material->alpha;
double beta = material->beta;
- double drosubg = 500.0 * context->gravity;
+ double drosub = 0.0f;
+ double p_est = context->pisos + 0.5f * ( phsDensity + drosub ) * context->gravity * ( radius - Spherical_RMin );
+ double rosubg = context->gravity * ( phsDensity + drosub ) * ( 1.0 - alpha * (nodeT-material->reftemp) + beta * p_est );
double press_norm = 0.0f;
float normal[3];
Normal* normal1;
@@ -550,7 +555,7 @@
mag_normal = sqrt( normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2] );
area *= 0.5f;
- press_norm = drosubg * ( Spherical_RMin - dhE );
+ press_norm = context->pisos + rosubg * ( Spherical_RMin - dhE );
(*force)[0] += factor4 * ( press_norm * area * normal[0] );
(*force)[1] += factor4 * ( press_norm * area * normal[1] );
(*force)[2] += factor4 * ( press_norm * area * normal[2] );
Modified: long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/InitialConditions.c 2008-10-15 19:11:30 UTC (rev 13088)
+++ long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/InitialConditions.c 2008-10-15 19:17:45 UTC (rev 13089)
@@ -46,24 +46,30 @@
Snac_Context* context = (Snac_Context*)_context;
- if( context->restartStep > 0 ) {
+ if( context->restartStep > 0 && (context->timeStep - context->restartStep == 1) ) {
FILE* fp;
char path[PATH_MAX];
Node_LocalIndex node_lI;
- Element_LocalIndex element_lI;
Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
sprintf( path, "%s/snac.isoForce.%d.restart", context->outputPath, context->rank );
fprintf(stderr,"%d: reading %s\n",context->rank,path);
- Journal_Firewall( ( ( fp = fopen( path, "r") ) != NULL ), context->snacError, "Failed to open %s!!\n", path );
+ Journal_Firewall( ( ( fp = fopen( path, "r") ) == NULL ), context->snacError, "Failed to open %s!!\n", path );
for( node_lI = 0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
Snac_Node* node = Snac_Node_At( context, node_lI );
fscanf( fp, "%le", &(node->residualFr) );
}
fclose( fp );
+
+ sprintf( path, "%s/pisos.restart", context->outputPath );
+ fprintf(stderr,"%d: reading %s\n",context->rank,path);
+ Journal_Firewall( ( ( fp = fopen( path, "r") ) == NULL ), context->snacError, "Failed to open %s!!\n", path );
+ fscanf( fp, "%le", &(context->pisos) );
+ fclose( fp );
+
}
}
Modified: long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/Register.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/Register.c 2008-10-15 19:11:30 UTC (rev 13088)
+++ long/3D/SNAC/trunk/Snac/plugins/winkler_G3_case3/Register.c 2008-10-15 19:17:45 UTC (rev 13089)
@@ -51,7 +51,6 @@
static double MIN[] = { -45.0f, -45.0f, 0.5f };
static double MAX[] = { 45.0f, 45.0f, 1.0f };
-
Index _SnacWinklerG3Force_Register( PluginsManager* pluginsMgr ) {
return PluginsManager_Submit( pluginsMgr,
SnacWinklerG3Force_Type,
@@ -75,9 +74,8 @@
name );
}
-
void _SnacWinklerG3Force_Construct( void* component, Stg_ComponentFactory* cf, void* data ) {
- Snac_Context* context;
+ Snac_Context* context;
Dictionary* meshDict;
int Spherical = 0;
Dictionary_Entry_Value* extensionsList;
More information about the cig-commits
mailing list