Improvment and bugfixing in CFiniteElementVec #105

wenqing · 2018-07-04T12:56:54Z

This PR improves CFiniteElementVec, and it contains

removed the negative value option of the Biot's constant.
in the assembly of the pressure coupling related RHS and matrix,
- replaced the element constant saturation with node saturations,
- changed the computation of pressure coupled matrix for the staggered scheme for H2M, and fixed a bug for the monolithic scheme of H2M
- removed the restriction of pore pressure to be positive in the total stress calculation.
fixed the excavation conditions for the time stepping excavation model.
replaced `int CFiniteElementVec::Flow_Type``with an enum variable.
made some variable declarations local.
removed some repeated or unused source code.

The changes in 2) improves the computation solutions of the benchmarks for the THM coupled processes. Therefore the reference results of the affected benchmarks have to be updated too.

wenqing · 2018-07-04T13:08:43Z

FEM/fem_ele_vec.cpp

 				{
 					if (smat->bishop_model < 0) // Without Bishop
 					{
 						for (int i = 0; i < nnodes; i++)
 						{
-							nodal_pore_p[i] = (biot < 0.0 && _nodal_p1[i] < 0.0) ?
-								 0.0 : LoadFactor * S_Water * _nodal_p1[i];
+							nodal_pore_p[i] = LoadFactor * _nodal_S[i] * _nodal_p1[i];


Element saturation to nodal saturation.

wenqing · 2018-07-04T13:10:09Z

FEM/fem_ele_vec.cpp

@@ -1721,7 +1683,7 @@ void CFiniteElementVec::GlobalAssembly_RHS()
 							{
 								const double p0 = h_pcs->GetNodeValue(nodes[i], idx_p1_ini);
 								const double Sat0 = LoadFactor * m_mmp->SaturationCapillaryPressureFunction(-p0);
-								nodal_pore_p[i] -= (p0 > 0.0) ? LoadFactor * Sat0 * p0 : 0.0;


removed the restriction of positive pore pressure in the total stress calculation for the partially saturated flow.

wenqing · 2018-07-04T13:11:13Z

FEM/fem_ele_std.cpp


 		for (i = 0; i < nnodesHQ; i++)
 		{
-			NodalVal2[i] = dm_pcs->GetNodeValue(nodes[i], Idx_dm1[0]);
-			NodalVal3[i] = dm_pcs->GetNodeValue(nodes[i], Idx_dm1[1]);
+			NodalVal2[i] = dm_pcs->GetNodeValue(nodes[i], Idx_dm0[0]);


Fixed a bug in the monolithic H2M.

wenqing · 2018-07-04T13:13:20Z

FEM/pcs_dm.cpp

@@ -2791,35 +2770,13 @@ void CRFProcessDeformation::PostExcavation()
 	{
 		//	 		  MXDumpGLS("rf_pcs.txt",1,eqs->b,eqs->x);  //abort();}
 		// 07.04.2010 WW
-		bool done;
-		ElementValue_DM* eleV_DM = NULL;
-		if ((fem_dm->Flow_Type == 0 || fem_dm->Flow_Type == 2 || fem_dm->Flow_Type == 1) && Neglect_H_ini == 2)


Removed the repeated code.

norihiro-w · 2018-07-10T23:17:19Z

FEM/fem_ele_vec.cpp

-					Sw += shapefct[i] * _nodal_S[i];
-			}
-			rho = (1. - porosity) * fabs(smat->Density()) + porosity * Sw * density_fluid;
+			const double density_fluid = (mfp_vector.empty()) ?  m_mfp->Density() : 0.0;


i am puzzled here. why should the fluid density be zero if mfp_vector is not empty?

Sorry, it is wrong. Removed the ternary operator.

norihiro-w · 2018-07-10T23:18:20Z

FEM/fem_ele_vec.cpp


-			if (Flow_Type == 2 || Flow_Type == 3)
+			double rho = (1. - porosity) * fabs(smat->Density())


since smat->Density() is positive (line 738), fabs is not needed.

norihiro-w · 2018-07-10T23:19:46Z

FEM/fem_ele_vec.cpp

-	return rho;
+
+    // If negative value is given in the .msp file, gravity by solid is skipped
+    if (smat->Density() > 0.0)


please correct indent

norihiro-w · 2018-07-10T23:25:08Z

FEM/fem_ele_vec.cpp

-			double fac = 1.0;
-			if (Flow_Type == 1)
-				fac = -1.0;
+			const double fac = (_flow_type == FiniteElement::RICHARDS_FLOW) ? -1. : 1.;


could you write a comment here explaining what is fac variable.

norihiro-w · 2018-07-10T23:32:19Z

FEM/rf_pcs.cpp

-			if ((GetCurveValue(ExcavCurve, 0, aktuelle_zeit, &valid) + ExcavBeginCoordinate)
-			        > (ele_center[ExcavDirection])
-			    && (ele_center[ExcavDirection] - ExcavBeginCoordinate) > -0.001)
+			const double expected_coordinate =


i see these lines several times. can you create a function for this and remove repeated codes?

Created a new function.

…ess results of TH_usM modelling

…CFiniteElementVec

in the the local assmebly.

…flow_type

to remove the duplicated pieces of code.

…lling

…ments

norihiro-w · 2018-08-02T21:11:03Z

FEM/fem_ele_vec.cpp

-				if ((GetCurveValue(pcs->ExcavCurve, 0, aktuelle_zeit, &valid) + pcs->ExcavBeginCoordinate)
-				        > (ele_center[pcs->ExcavDirection])
-				    && (ele_center[pcs->ExcavDirection] - pcs->ExcavBeginCoordinate) > -0.001)
+				const double expected_coordinate =


isn't it possible to use the new function isPointInExcavatedDomain() here?

Replace the lines with the new function.

norihiro-w · 2018-08-02T21:11:21Z

FEM/fem_ele_vec.cpp

@@ -1634,8 +1599,10 @@ void CFiniteElementVec::GlobalAssembly_RHS()
 				else if (pcs->ExcavMaterialGroup > -1)
 				{
 					double const* ele_center(elem->GetGravityCenter());
-					if ((GetCurveValue(pcs->ExcavCurve, 0, aktuelle_zeit, &valid) + pcs->ExcavBeginCoordinate)
-					    < (ele_center[pcs->ExcavDirection]))
+					const double expected_coordinate =


Just took a loot. Here the if-condition is different from that in the new function.

… code duplication

norihiro-w · 2018-08-07T07:23:42Z

looks good

This was referenced Jul 4, 2018

Update the results of the benchmarks for HM, THM and TH2M ufz/ogs5-benchmarks_ref#24

Merged

Modified the input files of THM and H2M_TEP ufz/ogs5-benchmarks#18

Merged

wenqing commented Jul 4, 2018

View reviewed changes

norihiro-w reviewed Jul 10, 2018

View reviewed changes

wenqing mentioned this pull request Jul 11, 2018

VTP + PVD output for particles if output option is set as PVD #106

Merged

wenqing mentioned this pull request Aug 1, 2018

[MPI/THM] Updated the results ufz/ogs5-benchmarks_ref#27

Merged

wenqing added 13 commits August 1, 2018 13:55

Removed the negative meaning of Biot constant, which improved the str…

1cebaf5

…ess results of TH_usM modelling

H2M: only compute one pressure coupled matrix for the staggered scheme

295a354

Changed the name and the calcualtion the wetting phase saturation in …

1322905

…CFiniteElementVec

H_usM: Repace the element constant saturation with node saturations

c95c0b1

in the the local assmebly.

CFiniteElementVec::Flow_Type: changed its type enum and its name to _…

6b1788c

…flow_type

CRFProcessDeformation: Use the existing member InitialNodeValueHpcs()

2eaf153

to remove the duplicated pieces of code.

Minor changes in the caculation of swelling models in CFiniteElementVec

00db87e

Little improvement in CFiniteElementVec::ComputeMatrix_RHS

2dc716d

Fixed the excavation conditions for the time stepping excavation mode…

f6d1054

…lling

Removed one compilation warning and corrected some indentations

5167f72

Fixed a bug in the assemble of RHS for the Newton scheme for H2M

c4744fd

[CFiniteElementVec] Comments and a correction according to Nori's com…

a4bc2d4

…ments

[CRFProcess] Added a member to check the excavation status of point

9fa57f0

wenqing force-pushed the improveLocalAssemblerM branch from 8992194 to c495882 Compare August 1, 2018 11:56

[CFiniteElementVec] Corrected spelling in comments

efe1ed9

wenqing force-pushed the improveLocalAssemblerM branch from c495882 to efe1ed9 Compare August 1, 2018 15:17

norihiro-w reviewed Aug 2, 2018

View reviewed changes

[CFiniteElementVec::LocalAssembly] Use isPointInExcavatedDomain avoid…

5943ad3

… code duplication

norihiro-w merged commit fe05142 into ufz:master Aug 7, 2018

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Improvment and bugfixing in CFiniteElementVec #105

Improvment and bugfixing in CFiniteElementVec #105

wenqing commented Jul 4, 2018 •

edited

Loading

wenqing Jul 4, 2018

wenqing Jul 4, 2018

wenqing Jul 4, 2018

wenqing Jul 4, 2018

norihiro-w Jul 10, 2018

wenqing Aug 1, 2018

norihiro-w Jul 10, 2018

wenqing Aug 1, 2018

norihiro-w Jul 10, 2018

wenqing Aug 1, 2018

norihiro-w Jul 10, 2018

wenqing Aug 1, 2018

norihiro-w Jul 10, 2018

wenqing Aug 1, 2018

norihiro-w Aug 2, 2018

wenqing Aug 3, 2018

norihiro-w Aug 2, 2018

wenqing Aug 3, 2018

norihiro-w commented Aug 7, 2018


		if (Flow_Type == 2 \|\| Flow_Type == 3)
		double rho = (1. - porosity) * fabs(smat->Density())

Improvment and bugfixing in CFiniteElementVec #105

Improvment and bugfixing in CFiniteElementVec #105

Conversation

wenqing commented Jul 4, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

norihiro-w commented Aug 7, 2018

wenqing commented Jul 4, 2018 •

edited

Loading