Sneddon manufactured solution. #1314
Conversation
isakhammer
requested review from
keileg,
IvarStefansson and
jwboth
as code owners
|
|
||
| bound_face_centers = g.face_centers[:, bound_faces] | ||
|
|
||
| for i in range(0, u_a.size): |
There was a problem hiding this comment.
Add some documentation.
| # Set the type of west and east boundaries to Dirichlet. North and south are | ||
| # Neumann by default. | ||
| bc = pp.BoundaryConditionVectorial(sd, bounds.all_bf, "dir") | ||
| frac_face = sd.tags["fracture_faces"] |
There was a problem hiding this comment.
I think there's a function for this. internal_bc_to_dirichlet or whatever.
There was a problem hiding this comment.
I also believe this to be the default behaviour - all Dirichlet on external and internal boundaries.
There was a problem hiding this comment.
Changed these lines to internal_to_dirichlet
tests/functional/test_sneddon_2d.py
Outdated
| """Test observed order of convergence. | ||
| """ | ||
| # We the order of L2 convergence on the fracture of displacement to be about 1.0 | ||
| assert 0.85 < actual_ooc["ooc_displacement"] |
There was a problem hiding this comment.
The test needs to be more precise.
| from porepy.geometry.distances import point_pointset | ||
|
|
||
|
|
||
| def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray: |
There was a problem hiding this comment.
Is a separate function needed for this, since it just wraps the other one?
| return cons * np.sqrt(1 - np.power(eta / a, 2)) | ||
|
|
||
|
|
||
| def transform(xc: np.ndarray, x: np.ndarray, alpha: float) -> np.ndarray: |
There was a problem hiding this comment.
Here and everywhere, would it be helpful to annotate the array sizes?
There was a problem hiding this comment.
I would rename this function, e.g. tranform_boundary_face_coordinates_bem
There was a problem hiding this comment.
The name may be improved, but the exact name may depend on context, i.e., Veljko's suggestion to move it.
| from porepy.geometry.distances import point_pointset | ||
|
|
||
|
|
||
| def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray: |
There was a problem hiding this comment.
This is a shallow wrapper which can be eliminated.
| return bem_centers | ||
|
|
||
|
|
||
| def analytical_displacements( |
There was a problem hiding this comment.
Should be moved to class representing the analytical solution.
| return cons * np.sqrt(1 - np.power(eta / a, 2)) | ||
|
|
||
|
|
||
| def transform(xc: np.ndarray, x: np.ndarray, alpha: float) -> np.ndarray: |
There was a problem hiding this comment.
Only used once in below method, can be eliminated and code put there.
In general, all these helper methods before the classes are created seem to be part of the analytical solution and can be moved there. Would also help reduce the signatures.
There was a problem hiding this comment.
Good observation. Please implement!
| u_bc = assign_bem(sd, h / 2, box_faces, self.theta, bem_centers, u_a, self.poi) | ||
| return u_bc | ||
|
|
||
| def exact_sol_fracture( |
There was a problem hiding this comment.
The first tuple element of return value serves currently no purpose anywhere.
|
|
||
| # Computing the numerical displacement jump along the fracture on the fracture | ||
| # cell centers. | ||
| u_n: pp.ad.Operator = nd_vec_to_normal @ self.displacement_jump(frac_sd) |
There was a problem hiding this comment.
type of variable u_n is typed explicitly, but then immediately changed in next line.
There was a problem hiding this comment.
Oh, dear. Rename below to u_n_value?
| u_n = u_n.value(self.equation_system) | ||
|
|
||
| # Checking convergence specifically on the fracture | ||
| u_a = self.exact_sol.exact_sol_fracture(self.mdg)[1] |
There was a problem hiding this comment.
Add truncation of around fracture tips to reach linear convergence (@IvarStefansson suggested 15% for example, but trial and error required).
There was a problem hiding this comment.
Please don't hard-code anything. I think we should allow the user to control this. Can be achieved by having a separate helper method and a parameter controlling the radius being removed around the tips.
tests/functional/test_sneddon_2d.py
Outdated
| } | ||
|
|
||
| # Convert angle to radians and compute fracture points | ||
| params["frac_pts"] = compute_frac_pts( |
There was a problem hiding this comment.
Since all the parameters required are given as model parameters, this computation can be moved to the custom geometry class of the setup (including the function used to compute the fracture points)
There was a problem hiding this comment.
Moved it, please check set_fractures.
|
|
||
| # Calculate and return the order of convergence for the displacement | ||
| order_dict = conv_analysis.order_of_convergence( | ||
| conv_analysis.run_analysis(), data_range=slice(None, None, None) |
There was a problem hiding this comment.
This is the default value for data_range.
| model_params=copy.deepcopy(params), | ||
| levels=2, | ||
| spatial_refinement_rate=2, | ||
| temporal_refinement_rate=1, |
There was a problem hiding this comment.
What's the point of this argument if the rate is 1? i.e. time step remains the same.
There was a problem hiding this comment.
True. The default is one, so we might as well remove it.
There was a problem hiding this comment.
Thanks for your efforts, @Yuriyzabegaev and @vlipovac. I have made some comments. Please implement everything you feel confident about, possibly consulting eachother and me.
| return bem_centers | ||
|
|
||
|
|
||
| def analytical_displacements( |
| return cons * np.sqrt(1 - np.power(eta / a, 2)) | ||
|
|
||
|
|
||
| def transform(xc: np.ndarray, x: np.ndarray, alpha: float) -> np.ndarray: |
There was a problem hiding this comment.
Good observation. Please implement!
| return cons * np.sqrt(1 - np.power(eta / a, 2)) | ||
|
|
||
|
|
||
| def transform(xc: np.ndarray, x: np.ndarray, alpha: float) -> np.ndarray: |
There was a problem hiding this comment.
The name may be improved, but the exact name may depend on context, i.e., Veljko's suggestion to move it.
|
|
||
| # Computing the numerical displacement jump along the fracture on the fracture | ||
| # cell centers. | ||
| u_n: pp.ad.Operator = nd_vec_to_normal @ self.displacement_jump(frac_sd) |
There was a problem hiding this comment.
Oh, dear. Rename below to u_n_value?
| u_n = u_n.value(self.equation_system) | ||
|
|
||
| # Checking convergence specifically on the fracture | ||
| u_a = self.exact_sol.exact_sol_fracture(self.mdg)[1] |
There was a problem hiding this comment.
Please don't hard-code anything. I think we should allow the user to control this. Can be achieved by having a separate helper method and a parameter controlling the radius being removed around the tips.
| model_params=copy.deepcopy(params), | ||
| levels=2, | ||
| spatial_refinement_rate=2, | ||
| temporal_refinement_rate=1, |
There was a problem hiding this comment.
True. The default is one, so we might as well remove it.
Proposed changes
This pull request improves upon the previously submitted PR (#1283). Runtime on my computer is approx 12 seconds.
This pull request includes reimplementing Sneddon's analytical solution as outlined in https://doi.org/10.1007/s10596-020-10002-5. Sneddon's solution provides an analytical method to compute the normal jump displacement for a pressurized fracture under linear elasticity conditions. The boundary element method is then integrated to compute semi-analytical boundary values, which are subsequently used to verify the MPSA numerical method for an open fracture. This process validates the accuracy and convergence of the proposed numerical method, making it well-suited as a test case.
Types of changes
What types of changes does this PR introduce to PorePy?
Put an
xin the boxes that apply.Checklist
Put an
xin the boxes that apply or explain briefly why the box is not relevant.pytestwas run with the--run-skippedflag.